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

Comparing:
branches/development/src/nonbonded/Electrostatic.cpp (file contents), Revision 1535 by gezelter, Fri Dec 31 18:31:56 2010 UTC vs.
branches/devel_omp/src/nonbonded/Electrostatic.cpp (file contents), Revision 1614 by mciznick, Tue Aug 23 20:55:51 2011 UTC

# Line 34 | Line 34
34   * work.  Good starting points are:
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 < * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
37 > * [2]  Fennell & Gezelter, J. Chem. Phys. 124 234104 (2006).          
38   * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39   * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
# Line 52 | Line 52 | namespace OpenMD {
52   namespace OpenMD {
53    
54    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
55 <                                  forceField_(NULL) {}
55 >                                  forceField_(NULL), info_(NULL),
56 >                                  haveCutoffRadius_(false),
57 >                                  haveDampingAlpha_(false),
58 >                                  haveDielectric_(false),
59 >                                  haveElectroSpline_(false)
60 >  {}
61    
62    void Electrostatic::initialize() {
63 +    
64 +    Globals* simParams_ = info_->getSimParams();
65  
59 –    Globals* simParams_;
60 –
66      summationMap_["HARD"]               = esm_HARD;
67      summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
68      summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
# Line 97 | Line 102 | namespace OpenMD {
102      screeningMethod_ = UNDAMPED;
103      dielectric_ = 1.0;
104      one_third_ = 1.0 / 3.0;
100 –    haveCutoffRadius_ = false;
101 –    haveDampingAlpha_ = false;
102 –    haveDielectric_ = false;  
103 –    haveElectroSpline_ = false;
105    
106      // check the summation method:
107      if (simParams_->haveElectrostaticSummationMethod()) {
# Line 113 | Line 114 | namespace OpenMD {
114        } else {
115          // throw error
116          sprintf( painCave.errMsg,
117 <                 "SimInfo error: Unknown electrostaticSummationMethod.\n"
117 >                 "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
118                   "\t(Input file specified %s .)\n"
119                   "\telectrostaticSummationMethod must be one of: \"none\",\n"
120                   "\t\"shifted_potential\", \"shifted_force\", or \n"
# Line 407 | Line 408 | namespace OpenMD {
408      return;
409    }
410    
411 <  void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
412 <                                                    RealType theRSW ) {
412 <    cutoffRadius_ = theECR;
411 >  void Electrostatic::setCutoffRadius( RealType rCut ) {
412 >    cutoffRadius_ = rCut;
413      rrf_ = cutoffRadius_;
414 –    rt_ = theRSW;
414      haveCutoffRadius_ = true;
415    }
416 +
417 +  void Electrostatic::setSwitchingRadius( RealType rSwitch ) {
418 +    rt_ = rSwitch;
419 +  }
420    void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
421      summationMethod_ = esm;
422    }
# Line 429 | Line 432 | namespace OpenMD {
432      haveDielectric_ = true;
433    }
434  
435 <  void Electrostatic::calcForce(InteractionData idat) {
435 >  void Electrostatic::initForce() {
436 >        if (!initialized_) initialize();
437 >  }
438  
439 +  void Electrostatic::calcForce(InteractionData &idat) {
440 +
441      // utility variables.  Should clean these up and use the Vector3d and
442      // Mat3x3d to replace as many as we can in future versions:
443  
# Line 443 | Line 450 | namespace OpenMD {
450      RealType ct_i, ct_j, ct_ij, a1;
451      RealType riji, ri, ri2, ri3, ri4;
452      RealType pref, vterm, epot, dudr;
453 +    RealType vpair(0.0);
454      RealType scale, sc2;
455      RealType pot_term, preVal, rfVal;
456      RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj;
457      RealType preSw, preSwSc;
458      RealType c1, c2, c3, c4;
459 <    RealType erfcVal, derfcVal;
459 >    RealType erfcVal(1.0), derfcVal(0.0);
460      RealType BigR;
461  
462      Vector3d Q_i, Q_j;
# Line 459 | Line 467 | namespace OpenMD {
467      Vector3d rhatdot2, rhatc4;
468      Vector3d dVdr;
469  
470 +    // variables for indirect (reaction field) interactions for excluded pairs:
471 +    RealType indirect_Pot(0.0);
472 +    RealType indirect_vpair(0.0);
473 +    Vector3d indirect_dVdr(V3Zero);
474 +    Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
475 +
476      pair<RealType, RealType> res;
477      
478      if (!initialized_) initialize();
479      
480 <    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1];
481 <    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2];
480 >    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
481 >    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
482      
483      // some variables we'll need independent of electrostatic type:
484  
485 <    riji = 1.0 / idat.rij;
486 <    Vector3d rhat = idat.d  * riji;
485 >    riji = 1.0 /  *(idat.rij) ;
486 >    Vector3d rhat =  *(idat.d)   * riji;
487  
488      // logicals
489  
# Line 483 | Line 497 | namespace OpenMD {
497      bool j_is_SplitDipole = data2.is_SplitDipole;
498      bool j_is_Quadrupole = data2.is_Quadrupole;
499      
500 <    if (i_is_Charge)
500 >    if (i_is_Charge) {
501        q_i = data1.charge;
502 +      if (idat.excluded) {
503 +        *(idat.skippedCharge2) += q_i;
504 +      }
505 +    }
506  
507      if (i_is_Dipole) {
508        mu_i = data1.dipole_moment;
509 <      uz_i = idat.eFrame1.getColumn(2);
509 >      uz_i = idat.eFrame1->getColumn(2);
510        
511        ct_i = dot(uz_i, rhat);
512  
# Line 504 | Line 522 | namespace OpenMD {
522        qyy_i = Q_i.y();
523        qzz_i = Q_i.z();
524        
525 <      ux_i = idat.eFrame1.getColumn(0);
526 <      uy_i = idat.eFrame1.getColumn(1);
527 <      uz_i = idat.eFrame1.getColumn(2);
525 >      ux_i = idat.eFrame1->getColumn(0);
526 >      uy_i = idat.eFrame1->getColumn(1);
527 >      uz_i = idat.eFrame1->getColumn(2);
528  
529        cx_i = dot(ux_i, rhat);
530        cy_i = dot(uy_i, rhat);
# Line 517 | Line 535 | namespace OpenMD {
535        duduz_i = V3Zero;
536      }
537  
538 <    if (j_is_Charge)
538 >    if (j_is_Charge) {
539        q_j = data2.charge;
540 +      if (idat.excluded) {
541 +        *(idat.skippedCharge1) += q_j;
542 +      }
543 +    }
544  
545 +
546      if (j_is_Dipole) {
547        mu_j = data2.dipole_moment;
548 <      uz_j = idat.eFrame2.getColumn(2);
548 >      uz_j = idat.eFrame2->getColumn(2);
549        
550        ct_j = dot(uz_j, rhat);
551  
# Line 538 | Line 561 | namespace OpenMD {
561        qyy_j = Q_j.y();
562        qzz_j = Q_j.z();
563        
564 <      ux_j = idat.eFrame2.getColumn(0);
565 <      uy_j = idat.eFrame2.getColumn(1);
566 <      uz_j = idat.eFrame2.getColumn(2);
564 >      ux_j = idat.eFrame2->getColumn(0);
565 >      uy_j = idat.eFrame2->getColumn(1);
566 >      uz_j = idat.eFrame2->getColumn(2);
567  
568        cx_j = dot(ux_j, rhat);
569        cy_j = dot(uy_j, rhat);
# Line 559 | Line 582 | namespace OpenMD {
582        if (j_is_Charge) {
583          if (screeningMethod_ == DAMPED) {
584            // assemble the damping variables
585 <          res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
585 >          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
586            erfcVal = res.first;
587            derfcVal = res.second;
588            c1 = erfcVal * riji;
# Line 569 | Line 592 | namespace OpenMD {
592            c2 = c1 * riji;
593          }
594  
595 <        preVal = idat.electroMult * pre11_ * q_i * q_j;
595 >        preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
596          
597          if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
598            vterm = preVal * (c1 - c1c_);
599 <          dudr  = -idat.sw * preVal * c2;
599 >          dudr  = - *(idat.sw)  * preVal * c2;
600  
601          } else if (summationMethod_ == esm_SHIFTED_FORCE)  {
602 <          vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - cutoffRadius_) );
603 <          dudr  = idat.sw * preVal * (c2c_ - c2);
602 >          vterm = preVal * ( c1 - c1c_ + c2c_*( *(idat.rij)  - cutoffRadius_) );
603 >          dudr  =  *(idat.sw)  * preVal * (c2c_ - c2);
604  
605          } else if (summationMethod_ == esm_REACTION_FIELD) {
606 <          rfVal = idat.electroMult * preRF_ * idat.rij * idat.rij;
606 >          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij);
607 >
608            vterm = preVal * ( riji + rfVal );            
609 <          dudr  = idat.sw * preVal * ( 2.0 * rfVal - riji ) * riji;
609 >          dudr  =  *(idat.sw)  * preVal * ( 2.0 * rfVal - riji ) * riji;
610 >          
611 >          // if this is an excluded pair, there are still indirect
612 >          // interactions via the reaction field we must worry about:
613  
614 <        } else {
615 <          vterm = preVal * riji * erfcVal;            
616 <
617 <          dudr  = - idat.sw * preVal * c2;
614 >          if (idat.excluded) {
615 >            indirect_vpair += preVal * rfVal;
616 >            indirect_Pot += *(idat.sw) * preVal * rfVal;
617 >            indirect_dVdr += *(idat.sw)  * preVal * 2.0 * rfVal  * riji * rhat;
618 >          }
619 >          
620 >        } else {
621  
622 +          vterm = preVal * riji * erfcVal;          
623 +          dudr  = -  *(idat.sw)  * preVal * c2;
624 +
625          }
593 –
594 –        idat.vpair += vterm;
595 –        epot += idat.sw * vterm;
626  
627 <        dVdr += dudr * rhat;      
627 >        vpair += vterm;
628 >        epot +=  *(idat.sw)  * vterm;
629 >        dVdr += dudr * rhat;                
630        }
631  
632        if (j_is_Dipole) {
633          // pref is used by all the possible methods
634 <        pref = idat.electroMult * pre12_ * q_i * mu_j;
635 <        preSw = idat.sw * pref;
634 >        pref =  *(idat.electroMult) * pre12_ * q_i * mu_j;
635 >        preSw =  *(idat.sw)  * pref;
636  
637          if (summationMethod_ == esm_REACTION_FIELD) {
638            ri2 = riji * riji;
639            ri3 = ri2 * riji;
640      
641 <          vterm = - pref * ct_j * ( ri2 - preRF2_ * idat.rij );
642 <          idat.vpair += vterm;
643 <          epot += idat.sw * vterm;
641 >          vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
642 >          vpair += vterm;
643 >          epot +=  *(idat.sw)  * vterm;
644  
645            dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
646 <          duduz_j += -preSw * rhat * (ri2 - preRF2_ * idat.rij);  
646 >          duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );  
647  
648 +          // Even if we excluded this pair from direct interactions,
649 +          // we still have the reaction-field-mediated charge-dipole
650 +          // interaction:
651 +
652 +          if (idat.excluded) {
653 +            indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij);
654 +            indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij);
655 +            indirect_dVdr += preSw * preRF2_ * uz_j;
656 +            indirect_duduz_j += preSw * rhat * preRF2_ *  *(idat.rij);
657 +          }
658 +                      
659          } else {
660            // determine the inverse r used if we have split dipoles
661            if (j_is_SplitDipole) {
662 <            BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
662 >            BigR = sqrt( *(idat.r2) + 0.25 * d_j * d_j);
663              ri = 1.0 / BigR;
664 <            scale = idat.rij * ri;
664 >            scale =  *(idat.rij)  * ri;
665            } else {
666              ri = riji;
667              scale = 1.0;
# Line 628 | Line 671 | namespace OpenMD {
671  
672            if (screeningMethod_ == DAMPED) {
673              // assemble the damping variables
674 <            res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
674 >            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
675              erfcVal = res.first;
676              derfcVal = res.second;
677              c1 = erfcVal * ri;
# Line 645 | Line 688 | namespace OpenMD {
688            // calculate the potential
689            pot_term =  scale * c2;
690            vterm = -pref * ct_j * pot_term;
691 <          idat.vpair += vterm;
692 <          epot += idat.sw * vterm;
691 >          vpair += vterm;
692 >          epot +=  *(idat.sw)  * vterm;
693              
694            // calculate derivatives for forces and torques
695  
# Line 661 | Line 704 | namespace OpenMD {
704          cx2 = cx_j * cx_j;
705          cy2 = cy_j * cy_j;
706          cz2 = cz_j * cz_j;
707 <        pref =  idat.electroMult * pre14_ * q_i * one_third_;
707 >        pref =   *(idat.electroMult) * pre14_ * q_i * one_third_;
708            
709          if (screeningMethod_ == DAMPED) {
710            // assemble the damping variables
711 <          res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
711 >          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
712            erfcVal = res.first;
713            derfcVal = res.second;
714            c1 = erfcVal * riji;
# Line 680 | Line 723 | namespace OpenMD {
723          }
724  
725          // precompute variables for convenience
726 <        preSw = idat.sw * pref;
726 >        preSw =  *(idat.sw)  * pref;
727          c2ri = c2 * riji;
728          c3ri = c3 * riji;
729 <        c4rij = c4 * idat.rij;
729 >        c4rij = c4 *  *(idat.rij) ;
730          rhatdot2 = 2.0 * rhat * c3;
731          rhatc4 = rhat * c4rij;
732  
# Line 692 | Line 735 | namespace OpenMD {
735                       qyy_j * (cy2*c3 - c2ri) +
736                       qzz_j * (cz2*c3 - c2ri) );
737          vterm = pref * pot_term;
738 <        idat.vpair += vterm;
739 <        epot += idat.sw * vterm;
738 >        vpair += vterm;
739 >        epot +=  *(idat.sw)  * vterm;
740                  
741          // calculate derivatives for the forces and torques
742  
# Line 711 | Line 754 | namespace OpenMD {
754  
755        if (j_is_Charge) {
756          // variables used by all the methods
757 <        pref = idat.electroMult * pre12_ * q_j * mu_i;
758 <        preSw = idat.sw * pref;
757 >        pref =  *(idat.electroMult) * pre12_ * q_j * mu_i;
758 >        preSw =  *(idat.sw)  * pref;
759  
760          if (summationMethod_ == esm_REACTION_FIELD) {
761  
762            ri2 = riji * riji;
763            ri3 = ri2 * riji;
764  
765 <          vterm = pref * ct_i * ( ri2 - preRF2_ * idat.rij );
766 <          idat.vpair += vterm;
767 <          epot += idat.sw * vterm;
765 >          vterm = pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
766 >          vpair += vterm;
767 >          epot +=  *(idat.sw)  * vterm;
768            
769            dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
770            
771 <          duduz_i += preSw * rhat * (ri2 - preRF2_ * idat.rij);
771 >          duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
772 >
773 >          // Even if we excluded this pair from direct interactions,
774 >          // we still have the reaction-field-mediated charge-dipole
775 >          // interaction:
776 >
777 >          if (idat.excluded) {
778 >            indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij);
779 >            indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij);
780 >            indirect_dVdr += -preSw * preRF2_ * uz_i;
781 >            indirect_duduz_i += -preSw * rhat * preRF2_ *  *(idat.rij);
782 >          }
783              
784          } else {
785            
786            // determine inverse r if we are using split dipoles
787            if (i_is_SplitDipole) {
788 <            BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
788 >            BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i);
789              ri = 1.0 / BigR;
790 <            scale = idat.rij * ri;
790 >            scale =  *(idat.rij)  * ri;
791            } else {
792              ri = riji;
793              scale = 1.0;
# Line 743 | Line 797 | namespace OpenMD {
797              
798            if (screeningMethod_ == DAMPED) {
799              // assemble the damping variables
800 <            res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
800 >            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
801              erfcVal = res.first;
802              derfcVal = res.second;
803              c1 = erfcVal * ri;
# Line 760 | Line 814 | namespace OpenMD {
814            // calculate the potential
815            pot_term = c2 * scale;
816            vterm = pref * ct_i * pot_term;
817 <          idat.vpair += vterm;
818 <          epot += idat.sw * vterm;
817 >          vpair += vterm;
818 >          epot +=  *(idat.sw)  * vterm;
819  
820            // calculate derivatives for the forces and torques
821            dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
# Line 773 | Line 827 | namespace OpenMD {
827          // variables used by all methods
828          ct_ij = dot(uz_i, uz_j);
829  
830 <        pref = idat.electroMult * pre22_ * mu_i * mu_j;
831 <        preSw = idat.sw * pref;
830 >        pref =  *(idat.electroMult) * pre22_ * mu_i * mu_j;
831 >        preSw =  *(idat.sw)  * pref;
832  
833          if (summationMethod_ == esm_REACTION_FIELD) {
834            ri2 = riji * riji;
# Line 783 | Line 837 | namespace OpenMD {
837  
838            vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
839                             preRF2_ * ct_ij );
840 <          idat.vpair += vterm;
841 <          epot += idat.sw * vterm;
840 >          vpair += vterm;
841 >          epot +=  *(idat.sw)  * vterm;
842              
843            a1 = 5.0 * ct_i * ct_j - ct_ij;
844              
# Line 793 | Line 847 | namespace OpenMD {
847            duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
848            duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
849  
850 +          if (idat.excluded) {
851 +            indirect_vpair +=  - pref * preRF2_ * ct_ij;
852 +            indirect_Pot +=    - preSw * preRF2_ * ct_ij;
853 +            indirect_duduz_i += -preSw * preRF2_ * uz_j;
854 +            indirect_duduz_j += -preSw * preRF2_ * uz_i;
855 +          }
856 +
857          } else {
858            
859            if (i_is_SplitDipole) {
860              if (j_is_SplitDipole) {
861 <              BigR = sqrt(idat.r2 + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
861 >              BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
862              } else {
863 <              BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
863 >              BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i);
864              }
865              ri = 1.0 / BigR;
866 <            scale = idat.rij * ri;
866 >            scale =  *(idat.rij)  * ri;
867            } else {
868              if (j_is_SplitDipole) {
869 <              BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
869 >              BigR = sqrt( *(idat.r2) + 0.25 * d_j * d_j);
870                ri = 1.0 / BigR;
871 <              scale = idat.rij * ri;
871 >              scale =  *(idat.rij)  * ri;
872              } else {
873                ri = riji;
874                scale = 1.0;
# Line 815 | Line 876 | namespace OpenMD {
876            }
877            if (screeningMethod_ == DAMPED) {
878              // assemble damping variables
879 <            res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
879 >            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
880              erfcVal = res.first;
881              derfcVal = res.second;
882              c1 = erfcVal * ri;
# Line 837 | Line 898 | namespace OpenMD {
898            preSwSc = preSw * scale;
899            c2ri = c2 * ri;
900            c3ri = c3 * ri;
901 <          c4rij = c4 * idat.rij;
901 >          c4rij = c4 *  *(idat.rij) ;
902  
903            // calculate the potential
904            pot_term = (ct_ij * c2ri - ctidotj * c3);
905            vterm = pref * pot_term;
906 <          idat.vpair += vterm;
907 <          epot += idat.sw * vterm;
906 >          vpair += vterm;
907 >          epot +=  *(idat.sw)  * vterm;
908  
909            // calculate derivatives for the forces and torques
910            dVdr += preSwSc * ( ctidotj * rhat * c4rij  -
# Line 862 | Line 923 | namespace OpenMD {
923          cy2 = cy_i * cy_i;
924          cz2 = cz_i * cz_i;
925  
926 <        pref = idat.electroMult * pre14_ * q_j * one_third_;
926 >        pref =  *(idat.electroMult) * pre14_ * q_j * one_third_;
927  
928          if (screeningMethod_ == DAMPED) {
929            // assemble the damping variables
930 <          res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
930 >          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
931            erfcVal = res.first;
932            derfcVal = res.second;
933            c1 = erfcVal * riji;
# Line 881 | Line 942 | namespace OpenMD {
942          }
943            
944          // precompute some variables for convenience
945 <        preSw = idat.sw * pref;
945 >        preSw =  *(idat.sw)  * pref;
946          c2ri = c2 * riji;
947          c3ri = c3 * riji;
948 <        c4rij = c4 * idat.rij;
948 >        c4rij = c4 *  *(idat.rij) ;
949          rhatdot2 = 2.0 * rhat * c3;
950          rhatc4 = rhat * c4rij;
951  
# Line 894 | Line 955 | namespace OpenMD {
955                       qzz_i * (cz2 * c3 - c2ri) );
956          
957          vterm = pref * pot_term;
958 <        idat.vpair += vterm;
959 <        epot += idat.sw * vterm;
958 >        vpair += vterm;
959 >        epot +=  *(idat.sw)  * vterm;
960  
961          // calculate the derivatives for the forces and torques
962  
# Line 909 | Line 970 | namespace OpenMD {
970        }
971      }
972  
912 –    idat.pot += epot;
913 –    idat.f1 += dVdr;
973  
974 <    if (i_is_Dipole || i_is_Quadrupole)
975 <      idat.t1 -= cross(uz_i, duduz_i);
976 <    if (i_is_Quadrupole) {
977 <      idat.t1 -= cross(ux_i, dudux_i);
978 <      idat.t1 -= cross(uy_i, duduy_i);
979 <    }
974 >    if (!idat.excluded) {
975 >      *(idat.vpair) += vpair;
976 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
977 >      *(idat.f1) += dVdr;
978 >      
979 >      if (i_is_Dipole || i_is_Quadrupole)
980 >        *(idat.t1) -= cross(uz_i, duduz_i);
981 >      if (i_is_Quadrupole) {
982 >        *(idat.t1) -= cross(ux_i, dudux_i);
983 >        *(idat.t1) -= cross(uy_i, duduy_i);
984 >      }
985 >      
986 >      if (j_is_Dipole || j_is_Quadrupole)
987 >        *(idat.t2) -= cross(uz_j, duduz_j);
988 >      if (j_is_Quadrupole) {
989 >        *(idat.t2) -= cross(uz_j, dudux_j);
990 >        *(idat.t2) -= cross(uz_j, duduy_j);
991 >      }
992  
993 <    if (j_is_Dipole || j_is_Quadrupole)
994 <      idat.t2 -= cross(uz_j, duduz_j);
995 <    if (j_is_Quadrupole) {
996 <      idat.t2 -= cross(uz_j, dudux_j);
997 <      idat.t2 -= cross(uz_j, duduy_j);
993 >    } else {
994 >
995 >      // only accumulate the forces and torques resulting from the
996 >      // indirect reaction field terms.
997 >      *(idat.vpair) += indirect_vpair;
998 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
999 >      *(idat.f1) += indirect_dVdr;
1000 >      
1001 >      if (i_is_Dipole)
1002 >        *(idat.t1) -= cross(uz_i, indirect_duduz_i);
1003 >      if (j_is_Dipole)
1004 >        *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1005      }
1006  
1007 +
1008      return;
1009    }  
1010  
1011 <  void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) {
1011 >  void Electrostatic::calcSkipCorrection(InteractionData &idat) {
1012  
1013      if (!initialized_) initialize();
1014      
1015 <    ElectrostaticAtomData data1 = ElectrostaticMap[skdat.atype1];
1016 <    ElectrostaticAtomData data2 = ElectrostaticMap[skdat.atype2];
1015 >    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
1016 >    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
1017      
1018      // logicals
1019  
# Line 945 | Line 1024 | namespace OpenMD {
1024      bool j_is_Dipole = data2.is_Dipole;
1025  
1026      RealType q_i, q_j;
1027 <    
1028 <    // The skippedCharge computation is needed by the real-space cutoff methods
1029 <    // (i.e. shifted force and shifted potential)
1030 <
1027 >
1028 >    // The skippedCharge computation is needed by the real-space
1029 >    // cutoff methods (i.e. shifted force and shifted potential)
1030 >      
1031      if (i_is_Charge) {
1032        q_i = data1.charge;
1033 <      skdat.skippedCharge2 += q_i;
1033 >      *(idat.skippedCharge2) += q_i;
1034      }
1035 <
1035 >    
1036      if (j_is_Charge) {
1037        q_j = data2.charge;
1038 <      skdat.skippedCharge1 += q_j;
1038 >      *(idat.skippedCharge1) += q_j;
1039      }
1040  
1041      // the rest of this function should only be necessary for reaction field.
1042  
1043      if (summationMethod_ == esm_REACTION_FIELD) {
1044        RealType riji, ri2, ri3;
1045 <      RealType q_i, mu_i, ct_i;
1046 <      RealType q_j, mu_j, ct_j;
1047 <      RealType preVal, rfVal, vterm, dudr, pref, myPot;
1045 >      RealType mu_i, ct_i;
1046 >      RealType mu_j, ct_j;
1047 >      RealType preVal, rfVal, vterm, dudr, pref, myPot(0.0);
1048        Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
1049  
1050        // some variables we'll need independent of electrostatic type:
1051        
1052 <      riji = 1.0 / skdat.rij;
1053 <      rhat = skdat.d  * riji;
1052 >      riji = 1.0 /  *(idat.rij) ;
1053 >      rhat =  *(idat.d)  * riji;
1054  
1055        if (i_is_Dipole) {
1056          mu_i = data1.dipole_moment;
1057 <        uz_i = skdat.eFrame1.getColumn(2);      
1057 >        uz_i = idat.eFrame1->getColumn(2);      
1058          ct_i = dot(uz_i, rhat);
1059          duduz_i = V3Zero;
1060        }
1061              
1062        if (j_is_Dipole) {
1063          mu_j = data2.dipole_moment;
1064 <        uz_j = skdat.eFrame2.getColumn(2);      
1064 >        uz_j = idat.eFrame2->getColumn(2);      
1065          ct_j = dot(uz_j, rhat);
1066          duduz_j = V3Zero;
1067        }
1068      
1069        if (i_is_Charge) {
1070          if (j_is_Charge) {
1071 <          preVal = skdat.electroMult * pre11_ * q_i * q_j;
1072 <          rfVal = preRF_ * skdat.rij * skdat.rij;
1071 >          preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
1072 >          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij) ;
1073            vterm = preVal * rfVal;
1074 <          myPot += skdat.sw * vterm;        
1075 <          dudr  = skdat.sw * preVal * 2.0 * rfVal * riji;        
1074 >          myPot +=  *(idat.sw)  * vterm;        
1075 >          dudr  =  *(idat.sw)  * preVal * 2.0 * rfVal * riji;        
1076            dVdr += dudr * rhat;
1077          }
1078          
1079          if (j_is_Dipole) {
1080            ri2 = riji * riji;
1081            ri3 = ri2 * riji;        
1082 <          pref = skdat.electroMult * pre12_ * q_i * mu_j;
1083 <          vterm = - pref * ct_j * ( ri2 - preRF2_ * skdat.rij );
1084 <          myPot += skdat.sw * vterm;        
1085 <          dVdr += -skdat.sw * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1086 <          duduz_j += -skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1082 >          pref =  *(idat.electroMult) * pre12_ * q_i * mu_j;
1083 >          vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
1084 >          myPot +=  *(idat.sw)  * vterm;        
1085 >          dVdr += - *(idat.sw)  * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1086 >          duduz_j += - *(idat.sw)  * pref * rhat * (ri2 - preRF2_ *  *(idat.rij) );
1087          }
1088        }
1089        if (i_is_Dipole) {
1090          if (j_is_Charge) {
1091            ri2 = riji * riji;
1092            ri3 = ri2 * riji;        
1093 <          pref = skdat.electroMult * pre12_ * q_j * mu_i;
1094 <          vterm = - pref * ct_i * ( ri2 - preRF2_ * skdat.rij );
1095 <          myPot += skdat.sw * vterm;        
1096 <          dVdr += skdat.sw * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);      
1097 <          duduz_i += skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1093 >          pref =  *(idat.electroMult) * pre12_ * q_j * mu_i;
1094 >          vterm = - pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
1095 >          myPot +=  *(idat.sw)  * vterm;        
1096 >          dVdr +=  *(idat.sw)  * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);      
1097 >          duduz_i +=  *(idat.sw)  * pref * rhat * (ri2 - preRF2_ *  *(idat.rij));
1098          }
1099        }
1100        
1101        // accumulate the forces and torques resulting from the self term
1102 <      skdat.pot += myPot;
1103 <      skdat.f1 += dVdr;
1102 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += myPot;
1103 >      *(idat.f1) += dVdr;
1104        
1105        if (i_is_Dipole)
1106 <        skdat.t1 -= cross(uz_i, duduz_i);
1106 >        *(idat.t1) -= cross(uz_i, duduz_i);
1107        if (j_is_Dipole)
1108 <        skdat.t2 -= cross(uz_j, duduz_j);
1108 >        *(idat.t2) -= cross(uz_j, duduz_j);
1109      }
1110    }
1111      
1112 <  void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) {
1112 >  void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1113      RealType mu1, preVal, chg1, self;
1114      
1115      if (!initialized_) initialize();
1116 <    
1117 <    ElectrostaticAtomData data = ElectrostaticMap[scdat.atype];
1116 >
1117 >    ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1118    
1119      // logicals
1041 –
1120      bool i_is_Charge = data.is_Charge;
1121      bool i_is_Dipole = data.is_Dipole;
1122  
# Line 1046 | Line 1124 | namespace OpenMD {
1124        if (i_is_Dipole) {
1125          mu1 = data.dipole_moment;          
1126          preVal = pre22_ * preRF2_ * mu1 * mu1;
1127 <        scdat.pot -= 0.5 * preVal;
1127 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal;
1128          
1129          // The self-correction term adds into the reaction field vector
1130 <        Vector3d uz_i = scdat.eFrame.getColumn(2);
1130 >        Vector3d uz_i = sdat.eFrame->getColumn(2);
1131          Vector3d ei = preVal * uz_i;
1132  
1133          // This looks very wrong.  A vector crossed with itself is zero.
1134 <        scdat.t -= cross(uz_i, ei);
1134 >        *(sdat.t) -= cross(uz_i, ei);
1135        }
1136      } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1137        if (i_is_Charge) {        
1138          chg1 = data.charge;
1139          if (screeningMethod_ == DAMPED) {
1140 <          self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1140 >          self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_;
1141          } else {        
1142 <          self = - 0.5 * rcuti_ * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1142 >          self = - 0.5 * rcuti_ * chg1 * (chg1 +  *(sdat.skippedCharge)) * pre11_;
1143          }
1144 <        scdat.pot += self;
1144 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1145        }
1146      }
1147    }
1148  
1149 <  RealType Electrostatic::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
1149 >  RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
1150      // This seems to work moderately well as a default.  There's no
1151      // inherent scale for 1/r interactions that we can standardize.
1152      // 12 angstroms seems to be a reasonably good guess for most

Diff Legend

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