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

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1584 by gezelter, Fri Jun 17 20:16:35 2011 UTC vs.
Revision 1668 by gezelter, Fri Jan 6 19:03:05 2012 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).                        
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include <stdio.h>
# Line 52 | Line 53 | namespace OpenMD {
53   namespace OpenMD {
54    
55    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
56 <                                  forceField_(NULL), info_(NULL) {}
56 >                                  forceField_(NULL), info_(NULL),
57 >                                  haveCutoffRadius_(false),
58 >                                  haveDampingAlpha_(false),
59 >                                  haveDielectric_(false),
60 >                                  haveElectroSpline_(false)
61 >  {}
62    
63    void Electrostatic::initialize() {
64 <
64 >    
65      Globals* simParams_ = info_->getSimParams();
66  
67      summationMap_["HARD"]               = esm_HARD;
68 +    summationMap_["NONE"]               = esm_HARD;
69      summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
70      summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
71      summationMap_["SHIFTED_FORCE"]      = esm_SHIFTED_FORCE;    
# Line 97 | Line 104 | namespace OpenMD {
104      screeningMethod_ = UNDAMPED;
105      dielectric_ = 1.0;
106      one_third_ = 1.0 / 3.0;
100 –    haveCutoffRadius_ = false;
101 –    haveDampingAlpha_ = false;
102 –    haveDielectric_ = false;  
103 –    haveElectroSpline_ = false;
107    
108      // check the summation method:
109      if (simParams_->haveElectrostaticSummationMethod()) {
# Line 115 | Line 118 | namespace OpenMD {
118          sprintf( painCave.errMsg,
119                   "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
120                   "\t(Input file specified %s .)\n"
121 <                 "\telectrostaticSummationMethod must be one of: \"none\",\n"
121 >                 "\telectrostaticSummationMethod must be one of: \"hard\",\n"
122                   "\t\"shifted_potential\", \"shifted_force\", or \n"
123                   "\t\"reaction_field\".\n", myMethod.c_str() );
124          painCave.isFatal = 1;
# Line 248 | Line 251 | namespace OpenMD {
251        preRF2_ = 2.0 * preRF_;
252      }
253      
254 <    RealType dx = cutoffRadius_ / RealType(np_ - 1);
254 >    // Add a 2 angstrom safety window to deal with cutoffGroups that
255 >    // have charged atoms longer than the cutoffRadius away from each
256 >    // other.  Splining may not be the best choice here.  Direct calls
257 >    // to erfc might be preferrable.
258 >
259 >    RealType dx = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
260      RealType rval;
261      vector<RealType> rvals;
262      vector<RealType> yvals;
# Line 445 | Line 453 | namespace OpenMD {
453      RealType ct_i, ct_j, ct_ij, a1;
454      RealType riji, ri, ri2, ri3, ri4;
455      RealType pref, vterm, epot, dudr;
456 +    RealType vpair(0.0);
457      RealType scale, sc2;
458      RealType pot_term, preVal, rfVal;
459      RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj;
460      RealType preSw, preSwSc;
461      RealType c1, c2, c3, c4;
462 <    RealType erfcVal, derfcVal;
462 >    RealType erfcVal(1.0), derfcVal(0.0);
463      RealType BigR;
464 +    RealType two(2.0), three(3.0);
465  
466      Vector3d Q_i, Q_j;
467      Vector3d ux_i, uy_i, uz_i;
# Line 460 | Line 470 | namespace OpenMD {
470      Vector3d dudux_j, duduy_j, duduz_j;
471      Vector3d rhatdot2, rhatc4;
472      Vector3d dVdr;
473 +
474 +    // variables for indirect (reaction field) interactions for excluded pairs:
475 +    RealType indirect_Pot(0.0);
476 +    RealType indirect_vpair(0.0);
477 +    Vector3d indirect_dVdr(V3Zero);
478 +    Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
479  
480      pair<RealType, RealType> res;
481      
# Line 485 | Line 501 | namespace OpenMD {
501      bool j_is_SplitDipole = data2.is_SplitDipole;
502      bool j_is_Quadrupole = data2.is_Quadrupole;
503      
504 <    if (i_is_Charge)
504 >    if (i_is_Charge) {
505        q_i = data1.charge;
506 +      if (idat.excluded) {
507 +        *(idat.skippedCharge2) += q_i;
508 +      }
509 +    }
510  
511      if (i_is_Dipole) {
512        mu_i = data1.dipole_moment;
# Line 519 | Line 539 | namespace OpenMD {
539        duduz_i = V3Zero;
540      }
541  
542 <    if (j_is_Charge)
542 >    if (j_is_Charge) {
543        q_j = data2.charge;
544 +      if (idat.excluded) {
545 +        *(idat.skippedCharge1) += q_j;
546 +      }
547 +    }
548  
549 +
550      if (j_is_Dipole) {
551        mu_j = data2.dipole_moment;
552        uz_j = idat.eFrame2->getColumn(2);
# Line 561 | Line 586 | namespace OpenMD {
586        if (j_is_Charge) {
587          if (screeningMethod_ == DAMPED) {
588            // assemble the damping variables
589 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
590 <          erfcVal = res.first;
591 <          derfcVal = res.second;
589 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
590 >          //erfcVal = res.first;
591 >          //derfcVal = res.second;
592 >
593 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
594 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
595 >
596            c1 = erfcVal * riji;
597            c2 = (-derfcVal + c1) * riji;
598          } else {
# Line 582 | Line 611 | namespace OpenMD {
611            dudr  =  *(idat.sw)  * preVal * (c2c_ - c2);
612  
613          } else if (summationMethod_ == esm_REACTION_FIELD) {
614 <          rfVal =  *(idat.electroMult) * preRF_ *  *(idat.rij)  *  *(idat.rij) ;
614 >          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij);
615 >
616            vterm = preVal * ( riji + rfVal );            
617            dudr  =  *(idat.sw)  * preVal * ( 2.0 * rfVal - riji ) * riji;
618 +          
619 +          // if this is an excluded pair, there are still indirect
620 +          // interactions via the reaction field we must worry about:
621  
622 +          if (idat.excluded) {
623 +            indirect_vpair += preVal * rfVal;
624 +            indirect_Pot += *(idat.sw) * preVal * rfVal;
625 +            indirect_dVdr += *(idat.sw)  * preVal * two * rfVal  * riji * rhat;
626 +          }
627 +          
628          } else {
590 –          vterm = preVal * riji * erfcVal;            
629  
630 +          vterm = preVal * riji * erfcVal;          
631            dudr  = -  *(idat.sw)  * preVal * c2;
632  
633          }
595 –
596 –        *(idat.vpair) += vterm;
597 –        epot +=  *(idat.sw)  * vterm;
634  
635 <        dVdr += dudr * rhat;      
635 >        vpair += vterm;
636 >        epot +=  *(idat.sw)  * vterm;
637 >        dVdr += dudr * rhat;                
638        }
639  
640        if (j_is_Dipole) {
# Line 609 | Line 647 | namespace OpenMD {
647            ri3 = ri2 * riji;
648      
649            vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
650 <          *(idat.vpair) += vterm;
650 >          vpair += vterm;
651            epot +=  *(idat.sw)  * vterm;
652  
653 <          dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
653 >          dVdr +=  -preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
654            duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );  
655  
656 +          // Even if we excluded this pair from direct interactions,
657 +          // we still have the reaction-field-mediated charge-dipole
658 +          // interaction:
659 +
660 +          if (idat.excluded) {
661 +            indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij);
662 +            indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij);
663 +            indirect_dVdr += preSw * preRF2_ * uz_j;
664 +            indirect_duduz_j += preSw * rhat * preRF2_ *  *(idat.rij);
665 +          }
666 +                      
667          } else {
668            // determine the inverse r used if we have split dipoles
669            if (j_is_SplitDipole) {
# Line 630 | Line 679 | namespace OpenMD {
679  
680            if (screeningMethod_ == DAMPED) {
681              // assemble the damping variables
682 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
683 <            erfcVal = res.first;
684 <            derfcVal = res.second;
682 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
683 >            //erfcVal = res.first;
684 >            //derfcVal = res.second;
685 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
686 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
687              c1 = erfcVal * ri;
688              c2 = (-derfcVal + c1) * ri;
689              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 647 | Line 698 | namespace OpenMD {
698            // calculate the potential
699            pot_term =  scale * c2;
700            vterm = -pref * ct_j * pot_term;
701 <          *(idat.vpair) += vterm;
701 >          vpair += vterm;
702            epot +=  *(idat.sw)  * vterm;
703              
704            // calculate derivatives for forces and torques
# Line 667 | Line 718 | namespace OpenMD {
718            
719          if (screeningMethod_ == DAMPED) {
720            // assemble the damping variables
721 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
722 <          erfcVal = res.first;
723 <          derfcVal = res.second;
721 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
722 >          //erfcVal = res.first;
723 >          //derfcVal = res.second;
724 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
725 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
726            c1 = erfcVal * riji;
727            c2 = (-derfcVal + c1) * riji;
728            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 686 | Line 739 | namespace OpenMD {
739          c2ri = c2 * riji;
740          c3ri = c3 * riji;
741          c4rij = c4 *  *(idat.rij) ;
742 <        rhatdot2 = 2.0 * rhat * c3;
742 >        rhatdot2 = two * rhat * c3;
743          rhatc4 = rhat * c4rij;
744  
745          // calculate the potential
# Line 694 | Line 747 | namespace OpenMD {
747                       qyy_j * (cy2*c3 - c2ri) +
748                       qzz_j * (cz2*c3 - c2ri) );
749          vterm = pref * pot_term;
750 <        *(idat.vpair) += vterm;
750 >        vpair += vterm;
751          epot +=  *(idat.sw)  * vterm;
752                  
753          // calculate derivatives for the forces and torques
754  
755 <        dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (2.0*cx_j*ux_j + rhat)*c3ri) +
756 <                           qyy_j* (cy2*rhatc4 - (2.0*cy_j*uy_j + rhat)*c3ri) +
757 <                           qzz_j* (cz2*rhatc4 - (2.0*cz_j*uz_j + rhat)*c3ri));
755 >        dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (two*cx_j*ux_j + rhat)*c3ri) +
756 >                           qyy_j* (cy2*rhatc4 - (two*cy_j*uy_j + rhat)*c3ri) +
757 >                           qzz_j* (cz2*rhatc4 - (two*cz_j*uz_j + rhat)*c3ri));
758                            
759          dudux_j += preSw * qxx_j * cx_j * rhatdot2;
760          duduy_j += preSw * qyy_j * cy_j * rhatdot2;
# Line 722 | Line 775 | namespace OpenMD {
775            ri3 = ri2 * riji;
776  
777            vterm = pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
778 <          *(idat.vpair) += vterm;
778 >          vpair += vterm;
779            epot +=  *(idat.sw)  * vterm;
780            
781 <          dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
781 >          dVdr += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_ * uz_i);
782            
783            duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
784 +
785 +          // Even if we excluded this pair from direct interactions,
786 +          // we still have the reaction-field-mediated charge-dipole
787 +          // interaction:
788 +
789 +          if (idat.excluded) {
790 +            indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij);
791 +            indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij);
792 +            indirect_dVdr += -preSw * preRF2_ * uz_i;
793 +            indirect_duduz_i += -preSw * rhat * preRF2_ *  *(idat.rij);
794 +          }
795              
796          } else {
797            
# Line 745 | Line 809 | namespace OpenMD {
809              
810            if (screeningMethod_ == DAMPED) {
811              // assemble the damping variables
812 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
813 <            erfcVal = res.first;
814 <            derfcVal = res.second;
812 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
813 >            //erfcVal = res.first;
814 >            //derfcVal = res.second;
815 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
816 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
817              c1 = erfcVal * ri;
818              c2 = (-derfcVal + c1) * ri;
819              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 762 | Line 828 | namespace OpenMD {
828            // calculate the potential
829            pot_term = c2 * scale;
830            vterm = pref * ct_i * pot_term;
831 <          *(idat.vpair) += vterm;
831 >          vpair += vterm;
832            epot +=  *(idat.sw)  * vterm;
833  
834            // calculate derivatives for the forces and torques
# Line 785 | Line 851 | namespace OpenMD {
851  
852            vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
853                             preRF2_ * ct_ij );
854 <          *(idat.vpair) += vterm;
854 >          vpair += vterm;
855            epot +=  *(idat.sw)  * vterm;
856              
857            a1 = 5.0 * ct_i * ct_j - ct_ij;
858              
859 <          dVdr += preSw * 3.0 * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
859 >          dVdr += preSw * three * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
860  
861 <          duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
862 <          duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
861 >          duduz_i += preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
862 >          duduz_j += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_*uz_i);
863  
864 +          if (idat.excluded) {
865 +            indirect_vpair +=  - pref * preRF2_ * ct_ij;
866 +            indirect_Pot +=    - preSw * preRF2_ * ct_ij;
867 +            indirect_duduz_i += -preSw * preRF2_ * uz_j;
868 +            indirect_duduz_j += -preSw * preRF2_ * uz_i;
869 +          }
870 +
871          } else {
872            
873            if (i_is_SplitDipole) {
# Line 817 | Line 890 | namespace OpenMD {
890            }
891            if (screeningMethod_ == DAMPED) {
892              // assemble damping variables
893 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
894 <            erfcVal = res.first;
895 <            derfcVal = res.second;
893 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
894 >            //erfcVal = res.first;
895 >            //derfcVal = res.second;
896 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
897 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
898              c1 = erfcVal * ri;
899              c2 = (-derfcVal + c1) * ri;
900              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 844 | Line 919 | namespace OpenMD {
919            // calculate the potential
920            pot_term = (ct_ij * c2ri - ctidotj * c3);
921            vterm = pref * pot_term;
922 <          *(idat.vpair) += vterm;
922 >          vpair += vterm;
923            epot +=  *(idat.sw)  * vterm;
924  
925            // calculate derivatives for the forces and torques
# Line 868 | Line 943 | namespace OpenMD {
943  
944          if (screeningMethod_ == DAMPED) {
945            // assemble the damping variables
946 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
947 <          erfcVal = res.first;
948 <          derfcVal = res.second;
946 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
947 >          //erfcVal = res.first;
948 >          //derfcVal = res.second;
949 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
950 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
951            c1 = erfcVal * riji;
952            c2 = (-derfcVal + c1) * riji;
953            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 887 | Line 964 | namespace OpenMD {
964          c2ri = c2 * riji;
965          c3ri = c3 * riji;
966          c4rij = c4 *  *(idat.rij) ;
967 <        rhatdot2 = 2.0 * rhat * c3;
967 >        rhatdot2 = two * rhat * c3;
968          rhatc4 = rhat * c4rij;
969  
970          // calculate the potential
# Line 896 | Line 973 | namespace OpenMD {
973                       qzz_i * (cz2 * c3 - c2ri) );
974          
975          vterm = pref * pot_term;
976 <        *(idat.vpair) += vterm;
976 >        vpair += vterm;
977          epot +=  *(idat.sw)  * vterm;
978  
979          // calculate the derivatives for the forces and torques
980  
981 <        dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (2.0*cx_i*ux_i + rhat)*c3ri) +
982 <                          qyy_i* (cy2*rhatc4 - (2.0*cy_i*uy_i + rhat)*c3ri) +
983 <                          qzz_i* (cz2*rhatc4 - (2.0*cz_i*uz_i + rhat)*c3ri));
981 >        dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (two*cx_i*ux_i + rhat)*c3ri) +
982 >                          qyy_i* (cy2*rhatc4 - (two*cy_i*uy_i + rhat)*c3ri) +
983 >                          qzz_i* (cz2*rhatc4 - (two*cz_i*uz_i + rhat)*c3ri));
984  
985          dudux_i += preSw * qxx_i * cx_i *  rhatdot2;
986          duduy_i += preSw * qyy_i * cy_i *  rhatdot2;
# Line 911 | Line 988 | namespace OpenMD {
988        }
989      }
990  
914 –    (*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
915 –    *(idat.f1) += dVdr;
991  
992 <    if (i_is_Dipole || i_is_Quadrupole)
993 <      *(idat.t1) -= cross(uz_i, duduz_i);
994 <    if (i_is_Quadrupole) {
995 <      *(idat.t1) -= cross(ux_i, dudux_i);
996 <      *(idat.t1) -= cross(uy_i, duduy_i);
997 <    }
998 <    
999 <    if (j_is_Dipole || j_is_Quadrupole)
1000 <      *(idat.t2) -= cross(uz_j, duduz_j);
1001 <    if (j_is_Quadrupole) {
1002 <      *(idat.t2) -= cross(uz_j, dudux_j);
1003 <      *(idat.t2) -= cross(uz_j, duduy_j);
1004 <    }
992 >    if (!idat.excluded) {
993 >      *(idat.vpair) += vpair;
994 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
995 >      *(idat.f1) += dVdr;
996 >      
997 >      if (i_is_Dipole || i_is_Quadrupole)
998 >        *(idat.t1) -= cross(uz_i, duduz_i);
999 >      if (i_is_Quadrupole) {
1000 >        *(idat.t1) -= cross(ux_i, dudux_i);
1001 >        *(idat.t1) -= cross(uy_i, duduy_i);
1002 >      }
1003 >      
1004 >      if (j_is_Dipole || j_is_Quadrupole)
1005 >        *(idat.t2) -= cross(uz_j, duduz_j);
1006 >      if (j_is_Quadrupole) {
1007 >        *(idat.t2) -= cross(uz_j, dudux_j);
1008 >        *(idat.t2) -= cross(uz_j, duduy_j);
1009 >      }
1010  
1011 <    return;
932 <  }  
1011 >    } else {
1012  
1013 <  void Electrostatic::calcSkipCorrection(InteractionData &idat) {
1013 >      // only accumulate the forces and torques resulting from the
1014 >      // indirect reaction field terms.
1015  
1016 <    if (!initialized_) initialize();
1017 <    
1018 <    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
939 <    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
940 <    
941 <    // logicals
942 <
943 <    bool i_is_Charge = data1.is_Charge;
944 <    bool i_is_Dipole = data1.is_Dipole;
945 <
946 <    bool j_is_Charge = data2.is_Charge;
947 <    bool j_is_Dipole = data2.is_Dipole;
948 <
949 <    RealType q_i, q_j;
950 <    
951 <    // The skippedCharge computation is needed by the real-space cutoff methods
952 <    // (i.e. shifted force and shifted potential)
953 <
954 <    if (i_is_Charge) {
955 <      q_i = data1.charge;
956 <      *(idat.skippedCharge2) += q_i;
957 <    }
958 <
959 <    if (j_is_Charge) {
960 <      q_j = data2.charge;
961 <      *(idat.skippedCharge1) += q_j;
962 <    }
963 <
964 <    // the rest of this function should only be necessary for reaction field.
965 <
966 <    if (summationMethod_ == esm_REACTION_FIELD) {
967 <      RealType riji, ri2, ri3;
968 <      RealType mu_i, ct_i;
969 <      RealType mu_j, ct_j;
970 <      RealType preVal, rfVal, vterm, dudr, pref, myPot(0.0);
971 <      Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
972 <
973 <      // some variables we'll need independent of electrostatic type:
1016 >      *(idat.vpair) += indirect_vpair;
1017 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
1018 >      *(idat.f1) += indirect_dVdr;
1019        
975 –      riji = 1.0 /  *(idat.rij) ;
976 –      rhat =  *(idat.d)  * riji;
977 –
978 –      if (i_is_Dipole) {
979 –        mu_i = data1.dipole_moment;
980 –        uz_i = idat.eFrame1->getColumn(2);      
981 –        ct_i = dot(uz_i, rhat);
982 –        duduz_i = V3Zero;
983 –      }
984 –            
985 –      if (j_is_Dipole) {
986 –        mu_j = data2.dipole_moment;
987 –        uz_j = idat.eFrame2->getColumn(2);      
988 –        ct_j = dot(uz_j, rhat);
989 –        duduz_j = V3Zero;
990 –      }
991 –    
992 –      if (i_is_Charge) {
993 –        if (j_is_Charge) {
994 –          preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
995 –          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij) ;
996 –          vterm = preVal * rfVal;
997 –          myPot +=  *(idat.sw)  * vterm;        
998 –          dudr  =  *(idat.sw)  * preVal * 2.0 * rfVal * riji;        
999 –          dVdr += dudr * rhat;
1000 –        }
1001 –        
1002 –        if (j_is_Dipole) {
1003 –          ri2 = riji * riji;
1004 –          ri3 = ri2 * riji;        
1005 –          pref =  *(idat.electroMult) * pre12_ * q_i * mu_j;
1006 –          vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
1007 –          myPot +=  *(idat.sw)  * vterm;        
1008 –          dVdr += - *(idat.sw)  * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1009 –          duduz_j += - *(idat.sw)  * pref * rhat * (ri2 - preRF2_ *  *(idat.rij) );
1010 –        }
1011 –      }
1012 –      if (i_is_Dipole) {
1013 –        if (j_is_Charge) {
1014 –          ri2 = riji * riji;
1015 –          ri3 = ri2 * riji;        
1016 –          pref =  *(idat.electroMult) * pre12_ * q_j * mu_i;
1017 –          vterm = - pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
1018 –          myPot +=  *(idat.sw)  * vterm;        
1019 –          dVdr +=  *(idat.sw)  * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);      
1020 –          duduz_i +=  *(idat.sw)  * pref * rhat * (ri2 - preRF2_ *  *(idat.rij));
1021 –        }
1022 –      }
1023 –      
1024 –      // accumulate the forces and torques resulting from the self term
1025 –      (*(idat.pot))[ELECTROSTATIC_FAMILY] += myPot;
1026 –      *(idat.f1) += dVdr;
1027 –      
1020        if (i_is_Dipole)
1021 <        *(idat.t1) -= cross(uz_i, duduz_i);
1021 >        *(idat.t1) -= cross(uz_i, indirect_duduz_i);
1022        if (j_is_Dipole)
1023 <        *(idat.t2) -= cross(uz_j, duduz_j);
1023 >        *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1024      }
1025 <  }
1025 >
1026 >
1027 >    return;
1028 >  }  
1029      
1030    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1031      RealType mu1, preVal, chg1, self;
1032      
1033      if (!initialized_) initialize();
1034 <    
1034 >
1035      ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1036    
1037      // logicals
1043 –
1038      bool i_is_Charge = data.is_Charge;
1039      bool i_is_Dipole = data.is_Dipole;
1040  
# Line 1048 | Line 1042 | namespace OpenMD {
1042        if (i_is_Dipole) {
1043          mu1 = data.dipole_moment;          
1044          preVal = pre22_ * preRF2_ * mu1 * mu1;
1045 <        sdat.pot[2] -= 0.5 * preVal;
1045 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal;
1046          
1047          // The self-correction term adds into the reaction field vector
1048          Vector3d uz_i = sdat.eFrame->getColumn(2);
# Line 1065 | Line 1059 | namespace OpenMD {
1059          } else {        
1060            self = - 0.5 * rcuti_ * chg1 * (chg1 +  *(sdat.skippedCharge)) * pre11_;
1061          }
1062 <        sdat.pot[ELECTROSTATIC_FAMILY] += self;
1062 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1063        }
1064      }
1065    }

Diff Legend

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