40 |
|
*/ |
41 |
|
|
42 |
|
#include "nonbonded/InteractionManager.hpp" |
43 |
– |
#include "UseTheForce/doForces_interface.h" |
43 |
|
|
44 |
|
namespace OpenMD { |
45 |
|
|
49 |
|
|
50 |
|
RealType InteractionManager::rCut_ = 0.0; |
51 |
|
RealType InteractionManager::rSwitch_ = 0.0; |
53 |
– |
RealType InteractionManager::skinThickness_ = 0.0; |
54 |
– |
RealType InteractionManager::listRadius_ = 0.0; |
52 |
|
CutoffMethod InteractionManager::cutoffMethod_ = SHIFTED_FORCE; |
53 |
|
SwitchingFunctionType InteractionManager::sft_ = cubic; |
54 |
|
RealType InteractionManager::vdwScale_[4] = {1.0, 0.0, 0.0, 0.0}; |
303 |
|
|
304 |
|
setupCutoffs(); |
305 |
|
setupSwitching(); |
309 |
– |
setupNeighborlists(); |
306 |
|
|
307 |
< |
int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; |
308 |
< |
int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; |
309 |
< |
notifyFortranCutoffs(&rCut_, &rSwitch_, &ljsp, &ljsf); |
314 |
< |
notifyFortranSkinThickness(&skinThickness_); |
307 |
> |
//int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; |
308 |
> |
//int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; |
309 |
> |
//notifyFortranCutoffs(&rCut_, &rSwitch_, &ljsp, &ljsf); |
310 |
|
|
311 |
|
initialized_ = true; |
312 |
|
} |
454 |
|
switcher_->setSwitch(rSwitch_, rCut_); |
455 |
|
} |
456 |
|
|
457 |
< |
/** |
463 |
< |
* setupNeighborlists |
464 |
< |
* |
465 |
< |
* If the skinThickness was explicitly set, use that value (but check it) |
466 |
< |
* If the skinThickness was not explicitly set: use 1.0 angstroms |
467 |
< |
*/ |
468 |
< |
void InteractionManager::setupNeighborlists() { |
469 |
< |
|
470 |
< |
Globals* simParams_ = info_->getSimParams(); |
471 |
< |
|
472 |
< |
if (simParams_->haveSkinThickness()) { |
473 |
< |
skinThickness_ = simParams_->getSkinThickness(); |
474 |
< |
} else { |
475 |
< |
skinThickness_ = 1.0; |
476 |
< |
sprintf(painCave.errMsg, |
477 |
< |
"InteractionManager::setupNeighborlists: No value was set for the skinThickness.\n" |
478 |
< |
"\tOpenMD will use a default value of %f Angstroms\n" |
479 |
< |
"\tfor this simulation\n", skinThickness_); |
480 |
< |
painCave.severity = OPENMD_INFO; |
481 |
< |
painCave.isFatal = 0; |
482 |
< |
simError(); |
483 |
< |
} |
484 |
< |
|
485 |
< |
listRadius_ = rCut_ + skinThickness_; |
486 |
< |
} |
487 |
< |
|
488 |
< |
|
489 |
< |
void InteractionManager::doPrePair(int *atid1, int *atid2, RealType *rij, RealType *rho_i_at_j, RealType *rho_j_at_i){ |
457 |
> |
void InteractionManager::doPrePair(InteractionData idat){ |
458 |
|
|
459 |
|
if (!initialized_) initialize(); |
460 |
< |
|
493 |
< |
DensityData ddat; |
494 |
< |
|
495 |
< |
ddat.atype1 = typeMap_[*atid1]; |
496 |
< |
ddat.atype2 = typeMap_[*atid2]; |
497 |
< |
ddat.rij = *rij; |
498 |
< |
ddat.rho_i_at_j = *rho_i_at_j; |
499 |
< |
ddat.rho_j_at_i = *rho_j_at_i; |
500 |
< |
|
501 |
< |
pair<AtomType*, AtomType*> key = make_pair(ddat.atype1, ddat.atype2); |
460 |
> |
|
461 |
|
set<NonBondedInteraction*>::iterator it; |
462 |
|
|
463 |
< |
for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ |
463 |
> |
for (it = interactions_[ *(idat.atypes) ].begin(); |
464 |
> |
it != interactions_[ *(idat.atypes) ].end(); ++it){ |
465 |
|
if ((*it)->getFamily() == METALLIC_FAMILY) { |
466 |
< |
dynamic_cast<MetallicInteraction*>(*it)->calcDensity(ddat); |
466 |
> |
dynamic_cast<MetallicInteraction*>(*it)->calcDensity(idat); |
467 |
|
} |
468 |
|
} |
469 |
|
|
470 |
|
return; |
471 |
|
} |
472 |
|
|
473 |
< |
void InteractionManager::doPreForce(int *atid, RealType *rho, RealType *frho, RealType *dfrhodrho){ |
473 |
> |
void InteractionManager::doPreForce(SelfData sdat){ |
474 |
|
|
475 |
|
if (!initialized_) initialize(); |
476 |
< |
|
477 |
< |
FunctionalData fdat; |
518 |
< |
|
519 |
< |
fdat.atype = typeMap_[*atid]; |
520 |
< |
fdat.rho = *rho; |
521 |
< |
fdat.frho = *frho; |
522 |
< |
fdat.dfrhodrho = *dfrhodrho; |
523 |
< |
|
524 |
< |
pair<AtomType*, AtomType*> key = make_pair(fdat.atype, fdat.atype); |
476 |
> |
|
477 |
> |
pair<AtomType*, AtomType*> key = make_pair(sdat.atype, sdat.atype); |
478 |
|
set<NonBondedInteraction*>::iterator it; |
479 |
|
|
480 |
|
for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ |
481 |
|
if ((*it)->getFamily() == METALLIC_FAMILY) { |
482 |
< |
dynamic_cast<MetallicInteraction*>(*it)->calcFunctional(fdat); |
482 |
> |
dynamic_cast<MetallicInteraction*>(*it)->calcFunctional(sdat); |
483 |
|
} |
484 |
|
} |
485 |
|
|
486 |
|
return; |
487 |
|
} |
488 |
|
|
489 |
< |
void InteractionManager::doPair(int *atid1, int *atid2, RealType *d, RealType *r, RealType *r2, RealType *sw, int *topoDist, RealType *A1, RealType *A2, RealType *eFrame1, RealType *eFrame2, RealType *vpair, RealType *pot, RealType *f1, RealType *t1, RealType *t2, RealType *rho1, RealType *rho2, RealType *dfrho1, RealType *dfrho2, RealType *fshift1, RealType *fshift2){ |
489 |
> |
void InteractionManager::doPair(InteractionData idat){ |
490 |
|
|
491 |
|
if (!initialized_) initialize(); |
492 |
< |
|
540 |
< |
InteractionData idat; |
541 |
< |
|
542 |
< |
idat.atype1 = typeMap_[*atid1]; |
543 |
< |
idat.atype2 = typeMap_[*atid2]; |
544 |
< |
idat.d = Vector3d(d); |
545 |
< |
idat.rij = *r; |
546 |
< |
idat.r2 = *r2; |
547 |
< |
idat.sw = *sw; |
548 |
< |
idat.vdwMult = vdwScale_[*topoDist]; |
549 |
< |
idat.electroMult = electrostaticScale_[*topoDist]; |
550 |
< |
for (int i = 0; i < 4; i++) { |
551 |
< |
idat.vpair[i] = vpair[i]; |
552 |
< |
idat.pot[i] = pot[i]; |
553 |
< |
} |
554 |
< |
idat.f1 = Vector3d(f1[0], f1[1], f1[2]); |
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); |
492 |
> |
|
493 |
|
set<NonBondedInteraction*>::iterator it; |
494 |
|
|
495 |
< |
for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) |
495 |
> |
for (it = interactions_[ *(idat.atypes) ].begin(); |
496 |
> |
it != interactions_[ *(idat.atypes) ].end(); ++it) |
497 |
|
(*it)->calcForce(idat); |
498 |
|
|
574 |
– |
for (int i = 0; i < 4; i++) { |
575 |
– |
vpair[i] = idat.vpair[i]; |
576 |
– |
pot[i] = idat.pot[i]; |
577 |
– |
} |
578 |
– |
|
579 |
– |
for (int i = 0; i < 3; i++) { |
580 |
– |
f1[i] = idat.f1[i]; |
581 |
– |
t1[i] = idat.t1[i]; |
582 |
– |
t2[i] = idat.t2[i]; |
583 |
– |
} |
584 |
– |
|
499 |
|
return; |
500 |
|
} |
501 |
|
|
502 |
< |
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){ |
502 |
> |
void InteractionManager::doSkipCorrection(InteractionData idat){ |
503 |
|
|
504 |
< |
if (!initialized_) initialize(); |
504 |
> |
if (!initialized_) initialize(); |
505 |
|
|
592 |
– |
SkipCorrectionData skdat; |
593 |
– |
|
594 |
– |
skdat.atype1 = typeMap_[*atid1]; |
595 |
– |
skdat.atype2 = typeMap_[*atid2]; |
596 |
– |
skdat.d = Vector3d(d); |
597 |
– |
skdat.rij = *r; |
598 |
– |
skdat.skippedCharge1 = *skippedCharge1; |
599 |
– |
skdat.skippedCharge2 = *skippedCharge2; |
600 |
– |
skdat.sw = *sw; |
601 |
– |
skdat.electroMult = *electroMult; |
602 |
– |
for (int i = 0; i < 4; i++) { |
603 |
– |
skdat.vpair[i] = vpair[i]; |
604 |
– |
skdat.pot[i] = pot[i]; |
605 |
– |
} |
606 |
– |
skdat.f1 = Vector3d(f1); |
607 |
– |
skdat.eFrame1 = Mat3x3d(eFrame1); |
608 |
– |
skdat.eFrame2 = Mat3x3d(eFrame2); |
609 |
– |
skdat.t1 = Vector3d(t1); |
610 |
– |
skdat.t2 = Vector3d(t2); |
611 |
– |
|
612 |
– |
pair<AtomType*, AtomType*> key = make_pair(skdat.atype1, skdat.atype2); |
506 |
|
set<NonBondedInteraction*>::iterator it; |
507 |
|
|
508 |
< |
for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ |
508 |
> |
for (it = interactions_[ *(idat.atypes) ].begin(); |
509 |
> |
it != interactions_[ *(idat.atypes) ].end(); ++it){ |
510 |
|
if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) { |
511 |
< |
dynamic_cast<ElectrostaticInteraction*>(*it)->calcSkipCorrection(skdat); |
511 |
> |
dynamic_cast<ElectrostaticInteraction*>(*it)->calcSkipCorrection(idat); |
512 |
|
} |
513 |
|
} |
514 |
|
|
621 |
– |
for (int i = 0; i < 4; i++) { |
622 |
– |
vpair[i] = skdat.vpair[i]; |
623 |
– |
pot[i] = skdat.pot[i]; |
624 |
– |
} |
625 |
– |
|
626 |
– |
for (int i = 0; i < 3; i++) { |
627 |
– |
f1[i] = skdat.f1[i]; |
628 |
– |
t1[i] = skdat.t1[i]; |
629 |
– |
t2[i] = skdat.t2[i]; |
630 |
– |
} |
631 |
– |
|
515 |
|
return; |
516 |
|
} |
517 |
|
|
518 |
< |
void InteractionManager::doSelfCorrection(int *atid, RealType *eFrame, RealType *skippedCharge, RealType *pot, RealType *t){ |
518 |
> |
void InteractionManager::doSelfCorrection(SelfData sdat){ |
519 |
|
|
520 |
|
if (!initialized_) initialize(); |
521 |
|
|
522 |
< |
SelfCorrectionData scdat; |
640 |
< |
|
641 |
< |
scdat.atype = typeMap_[*atid]; |
642 |
< |
scdat.eFrame = Mat3x3d(eFrame); |
643 |
< |
scdat.skippedCharge = *skippedCharge; |
644 |
< |
for (int i = 0; i < 4; i++){ |
645 |
< |
scdat.pot[i] = pot[i]; |
646 |
< |
} |
647 |
< |
scdat.t = Vector3d(t); |
648 |
< |
|
649 |
< |
pair<AtomType*, AtomType*> key = make_pair(scdat.atype, scdat.atype); |
522 |
> |
pair<AtomType*, AtomType*> key = make_pair(sdat.atype, sdat.atype); |
523 |
|
set<NonBondedInteraction*>::iterator it; |
524 |
|
|
525 |
|
for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ |
526 |
|
if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) { |
527 |
< |
dynamic_cast<ElectrostaticInteraction*>(*it)->calcSelfCorrection(scdat); |
527 |
> |
dynamic_cast<ElectrostaticInteraction*>(*it)->calcSelfCorrection(sdat); |
528 |
|
} |
529 |
|
} |
530 |
|
|
658 |
– |
t[0] = scdat.t.x(); |
659 |
– |
t[1] = scdat.t.y(); |
660 |
– |
t[2] = scdat.t.z(); |
661 |
– |
|
531 |
|
return; |
532 |
|
} |
533 |
|
|
665 |
– |
|
534 |
|
RealType InteractionManager::getSuggestedCutoffRadius(int *atid) { |
535 |
|
if (!initialized_) initialize(); |
536 |
|
|
541 |
|
RealType cutoff = 0.0; |
542 |
|
|
543 |
|
for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) |
544 |
< |
cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(atype, atype)); |
544 |
> |
cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key)); |
545 |
|
return cutoff; |
546 |
|
} |
547 |
|
|
553 |
|
RealType cutoff = 0.0; |
554 |
|
|
555 |
|
for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) |
556 |
< |
cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(atype, atype)); |
556 |
> |
cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key)); |
557 |
|
return cutoff; |
558 |
|
} |
691 |
– |
|
692 |
– |
void InteractionManager::getSwitch(RealType *r2, RealType *sw, RealType *dswdr, RealType *r, |
693 |
– |
int *in_switching_region) { |
694 |
– |
bool isr = switcher_->getSwitch( *r2 , *sw, *dswdr, *r); |
695 |
– |
*in_switching_region = (int)isr; |
696 |
– |
} |
697 |
– |
|
559 |
|
} //end namespace OpenMD |
699 |
– |
|
700 |
– |
extern "C" { |
701 |
– |
|
702 |
– |
#define fortranDoPrePair FC_FUNC(do_prepair, DO_PREPAIR) |
703 |
– |
#define fortranDoPreForce FC_FUNC(do_preforce, DO_PREFORCE) |
704 |
– |
#define fortranDoPair FC_FUNC(do_pair, DO_PAIR) |
705 |
– |
#define fortranDoSkipCorrection FC_FUNC(do_skip_correction, DO_SKIP_CORRECTION) |
706 |
– |
#define fortranDoSelfCorrection FC_FUNC(do_self_correction, DO_SELF_CORRECTION) |
707 |
– |
#define fortranGetCutoff FC_FUNC(get_cutoff, GET_CUTOFF) |
708 |
– |
#define fortranGetSwitch FC_FUNC(get_switch, GET_SWITCH) |
709 |
– |
|
710 |
– |
void fortranDoPrePair(int *atid1, int *atid2, RealType *rij, |
711 |
– |
RealType *rho_i_at_j, RealType *rho_j_at_i) { |
712 |
– |
|
713 |
– |
return OpenMD::InteractionManager::Instance()->doPrePair(atid1, atid2, rij, |
714 |
– |
rho_i_at_j, |
715 |
– |
rho_j_at_i); |
716 |
– |
} |
717 |
– |
void fortranDoPreForce(int *atid, RealType *rho, RealType *frho, |
718 |
– |
RealType *dfrhodrho) { |
719 |
– |
|
720 |
– |
return OpenMD::InteractionManager::Instance()->doPreForce(atid, rho, frho, |
721 |
– |
dfrhodrho); |
722 |
– |
} |
723 |
– |
|
724 |
– |
void fortranDoPair(int *atid1, int *atid2, RealType *d, RealType *r, RealType *r2, |
725 |
– |
RealType *sw, int *topoDist, RealType *A1, RealType *A2, |
726 |
– |
RealType *eFrame1, RealType *eFrame2, |
727 |
– |
RealType *vpair, RealType *pot, |
728 |
– |
RealType *f1, RealType *t1, RealType *t2, |
729 |
– |
RealType *rho1, RealType *rho2, RealType *dfrho1, RealType *dfrho2, |
730 |
– |
RealType *fshift1, RealType *fshift2){ |
731 |
– |
|
732 |
– |
return OpenMD::InteractionManager::Instance()->doPair(atid1, atid2, d, r, r2, |
733 |
– |
sw, topoDist, A1, A2, |
734 |
– |
eFrame1, eFrame2, |
735 |
– |
vpair, pot, |
736 |
– |
f1, t1, t2, |
737 |
– |
rho1, rho2, |
738 |
– |
dfrho1, dfrho2, |
739 |
– |
fshift1, fshift2); |
740 |
– |
} |
741 |
– |
|
742 |
– |
void fortranDoSkipCorrection(int *atid1, int *atid2, RealType *d, |
743 |
– |
RealType *r, RealType *skippedCharge1, |
744 |
– |
RealType *skippedCharge2, RealType *sw, |
745 |
– |
RealType *electroMult, RealType *pot, |
746 |
– |
RealType *vpair, RealType *f1, |
747 |
– |
RealType *eFrame1, RealType *eFrame2, |
748 |
– |
RealType *t1, RealType *t2){ |
749 |
– |
|
750 |
– |
return OpenMD::InteractionManager::Instance()->doSkipCorrection(atid1, |
751 |
– |
atid2, d, |
752 |
– |
r, |
753 |
– |
skippedCharge1, |
754 |
– |
skippedCharge2, |
755 |
– |
sw, electroMult, pot, |
756 |
– |
vpair, f1, eFrame1, |
757 |
– |
eFrame2, t1, t2); |
758 |
– |
} |
759 |
– |
|
760 |
– |
void fortranDoSelfCorrection(int *atid, RealType *eFrame, RealType *skippedCharge, |
761 |
– |
RealType *pot, RealType *t) { |
762 |
– |
|
763 |
– |
return OpenMD::InteractionManager::Instance()->doSelfCorrection(atid, |
764 |
– |
eFrame, |
765 |
– |
skippedCharge, |
766 |
– |
pot, t); |
767 |
– |
} |
768 |
– |
RealType fortranGetCutoff(int *atid) { |
769 |
– |
return OpenMD::InteractionManager::Instance()->getSuggestedCutoffRadius(atid); |
770 |
– |
} |
771 |
– |
|
772 |
– |
void fortranGetSwitch(RealType *r2, RealType *sw, RealType *dswdr, RealType *r, |
773 |
– |
int *in_switching_region) { |
774 |
– |
return OpenMD::InteractionManager::Instance()->getSwitch(r2, sw, dswdr, r, |
775 |
– |
in_switching_region); |
776 |
– |
} |
777 |
– |
} |