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 1554 by gezelter, Sat Apr 30 02:54:02 2011 UTC vs.
Revision 1665 by gezelter, Tue Nov 22 20:38:56 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).                        
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) {}
56 >                                  forceField_(NULL), info_(NULL),
57 >                                  haveCutoffRadius_(false),
58 >                                  haveDampingAlpha_(false),
59 >                                  haveDielectric_(false),
60 >                                  haveElectroSpline_(false)
61 >  {}
62    
63    void Electrostatic::initialize() {
64 +    
65 +    Globals* simParams_ = info_->getSimParams();
66  
59 –    Globals* simParams_;
60 –
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 407 | Line 415 | namespace OpenMD {
415      return;
416    }
417    
418 <  void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
419 <                                                    RealType theRSW ) {
412 <    cutoffRadius_ = theECR;
418 >  void Electrostatic::setCutoffRadius( RealType rCut ) {
419 >    cutoffRadius_ = rCut;
420      rrf_ = cutoffRadius_;
414 –    rt_ = theRSW;
421      haveCutoffRadius_ = true;
422    }
423 +
424 +  void Electrostatic::setSwitchingRadius( RealType rSwitch ) {
425 +    rt_ = rSwitch;
426 +  }
427    void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
428      summationMethod_ = esm;
429    }
# Line 443 | 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  
465      Vector3d Q_i, Q_j;
# Line 459 | Line 470 | namespace OpenMD {
470      Vector3d rhatdot2, rhatc4;
471      Vector3d dVdr;
472  
473 +    // variables for indirect (reaction field) interactions for excluded pairs:
474 +    RealType indirect_Pot(0.0);
475 +    RealType indirect_vpair(0.0);
476 +    Vector3d indirect_dVdr(V3Zero);
477 +    Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
478 +
479      pair<RealType, RealType> res;
480      
481      if (!initialized_) initialize();
482      
483 <    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes->first];
484 <    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes->second];
483 >    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
484 >    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
485      
486      // some variables we'll need independent of electrostatic type:
487  
# Line 483 | Line 500 | namespace OpenMD {
500      bool j_is_SplitDipole = data2.is_SplitDipole;
501      bool j_is_Quadrupole = data2.is_Quadrupole;
502      
503 <    if (i_is_Charge)
503 >    if (i_is_Charge) {
504        q_i = data1.charge;
505 +      if (idat.excluded) {
506 +        *(idat.skippedCharge2) += q_i;
507 +      }
508 +    }
509  
510      if (i_is_Dipole) {
511        mu_i = data1.dipole_moment;
# Line 517 | Line 538 | namespace OpenMD {
538        duduz_i = V3Zero;
539      }
540  
541 <    if (j_is_Charge)
541 >    if (j_is_Charge) {
542        q_j = data2.charge;
543 +      if (idat.excluded) {
544 +        *(idat.skippedCharge1) += q_j;
545 +      }
546 +    }
547  
548 +
549      if (j_is_Dipole) {
550        mu_j = data2.dipole_moment;
551        uz_j = idat.eFrame2->getColumn(2);
# Line 559 | Line 585 | namespace OpenMD {
585        if (j_is_Charge) {
586          if (screeningMethod_ == DAMPED) {
587            // assemble the damping variables
588 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
589 <          erfcVal = res.first;
590 <          derfcVal = res.second;
588 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
589 >          //erfcVal = res.first;
590 >          //derfcVal = res.second;
591 >
592 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
593 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
594 >
595            c1 = erfcVal * riji;
596            c2 = (-derfcVal + c1) * riji;
597          } else {
# Line 580 | Line 610 | namespace OpenMD {
610            dudr  =  *(idat.sw)  * preVal * (c2c_ - c2);
611  
612          } else if (summationMethod_ == esm_REACTION_FIELD) {
613 <          rfVal =  *(idat.electroMult) * preRF_ *  *(idat.rij)  *  *(idat.rij) ;
613 >          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij);
614 >
615            vterm = preVal * ( riji + rfVal );            
616            dudr  =  *(idat.sw)  * preVal * ( 2.0 * rfVal - riji ) * riji;
617 +          
618 +          // if this is an excluded pair, there are still indirect
619 +          // interactions via the reaction field we must worry about:
620  
621 <        } else {
622 <          vterm = preVal * riji * erfcVal;            
621 >          if (idat.excluded) {
622 >            indirect_vpair += preVal * rfVal;
623 >            indirect_Pot += *(idat.sw) * preVal * rfVal;
624 >            indirect_dVdr += *(idat.sw)  * preVal * 2.0 * rfVal  * riji * rhat;
625 >          }
626 >          
627 >        } else {
628  
629 +          vterm = preVal * riji * erfcVal;          
630            dudr  = -  *(idat.sw)  * preVal * c2;
631  
632          }
593 –
594 –        *(idat.vpair) += vterm;
595 –        epot +=  *(idat.sw)  * vterm;
633  
634 <        dVdr += dudr * rhat;      
634 >        vpair += vterm;
635 >        epot +=  *(idat.sw)  * vterm;
636 >        dVdr += dudr * rhat;                
637        }
638  
639        if (j_is_Dipole) {
# Line 607 | Line 646 | namespace OpenMD {
646            ri3 = ri2 * riji;
647      
648            vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
649 <          *(idat.vpair) += vterm;
649 >          vpair += vterm;
650            epot +=  *(idat.sw)  * vterm;
651  
652            dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
653            duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );  
654  
655 +          // Even if we excluded this pair from direct interactions,
656 +          // we still have the reaction-field-mediated charge-dipole
657 +          // interaction:
658 +
659 +          if (idat.excluded) {
660 +            indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij);
661 +            indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij);
662 +            indirect_dVdr += preSw * preRF2_ * uz_j;
663 +            indirect_duduz_j += preSw * rhat * preRF2_ *  *(idat.rij);
664 +          }
665 +                      
666          } else {
667            // determine the inverse r used if we have split dipoles
668            if (j_is_SplitDipole) {
# Line 628 | Line 678 | namespace OpenMD {
678  
679            if (screeningMethod_ == DAMPED) {
680              // assemble the damping variables
681 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
682 <            erfcVal = res.first;
683 <            derfcVal = res.second;
681 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
682 >            //erfcVal = res.first;
683 >            //derfcVal = res.second;
684 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
685 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
686              c1 = erfcVal * ri;
687              c2 = (-derfcVal + c1) * ri;
688              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 645 | Line 697 | namespace OpenMD {
697            // calculate the potential
698            pot_term =  scale * c2;
699            vterm = -pref * ct_j * pot_term;
700 <          *(idat.vpair) += vterm;
700 >          vpair += vterm;
701            epot +=  *(idat.sw)  * vterm;
702              
703            // calculate derivatives for forces and torques
# Line 665 | Line 717 | namespace OpenMD {
717            
718          if (screeningMethod_ == DAMPED) {
719            // assemble the damping variables
720 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
721 <          erfcVal = res.first;
722 <          derfcVal = res.second;
720 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
721 >          //erfcVal = res.first;
722 >          //derfcVal = res.second;
723 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
724 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
725            c1 = erfcVal * riji;
726            c2 = (-derfcVal + c1) * riji;
727            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 692 | Line 746 | namespace OpenMD {
746                       qyy_j * (cy2*c3 - c2ri) +
747                       qzz_j * (cz2*c3 - c2ri) );
748          vterm = pref * pot_term;
749 <        *(idat.vpair) += vterm;
749 >        vpair += vterm;
750          epot +=  *(idat.sw)  * vterm;
751                  
752          // calculate derivatives for the forces and torques
# Line 720 | Line 774 | namespace OpenMD {
774            ri3 = ri2 * riji;
775  
776            vterm = pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
777 <          *(idat.vpair) += vterm;
777 >          vpair += vterm;
778            epot +=  *(idat.sw)  * vterm;
779            
780            dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
781            
782            duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
783 +
784 +          // Even if we excluded this pair from direct interactions,
785 +          // we still have the reaction-field-mediated charge-dipole
786 +          // interaction:
787 +
788 +          if (idat.excluded) {
789 +            indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij);
790 +            indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij);
791 +            indirect_dVdr += -preSw * preRF2_ * uz_i;
792 +            indirect_duduz_i += -preSw * rhat * preRF2_ *  *(idat.rij);
793 +          }
794              
795          } else {
796            
# Line 743 | Line 808 | namespace OpenMD {
808              
809            if (screeningMethod_ == DAMPED) {
810              // assemble the damping variables
811 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
812 <            erfcVal = res.first;
813 <            derfcVal = res.second;
811 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
812 >            //erfcVal = res.first;
813 >            //derfcVal = res.second;
814 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
815 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
816              c1 = erfcVal * ri;
817              c2 = (-derfcVal + c1) * ri;
818              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 760 | Line 827 | namespace OpenMD {
827            // calculate the potential
828            pot_term = c2 * scale;
829            vterm = pref * ct_i * pot_term;
830 <          *(idat.vpair) += vterm;
830 >          vpair += vterm;
831            epot +=  *(idat.sw)  * vterm;
832  
833            // calculate derivatives for the forces and torques
# Line 783 | Line 850 | namespace OpenMD {
850  
851            vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
852                             preRF2_ * ct_ij );
853 <          *(idat.vpair) += vterm;
853 >          vpair += vterm;
854            epot +=  *(idat.sw)  * vterm;
855              
856            a1 = 5.0 * ct_i * ct_j - ct_ij;
# Line 793 | Line 860 | namespace OpenMD {
860            duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
861            duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
862  
863 +          if (idat.excluded) {
864 +            indirect_vpair +=  - pref * preRF2_ * ct_ij;
865 +            indirect_Pot +=    - preSw * preRF2_ * ct_ij;
866 +            indirect_duduz_i += -preSw * preRF2_ * uz_j;
867 +            indirect_duduz_j += -preSw * preRF2_ * uz_i;
868 +          }
869 +
870          } else {
871            
872            if (i_is_SplitDipole) {
# Line 815 | Line 889 | namespace OpenMD {
889            }
890            if (screeningMethod_ == DAMPED) {
891              // assemble damping variables
892 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
893 <            erfcVal = res.first;
894 <            derfcVal = res.second;
892 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
893 >            //erfcVal = res.first;
894 >            //derfcVal = res.second;
895 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
896 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
897              c1 = erfcVal * ri;
898              c2 = (-derfcVal + c1) * ri;
899              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 842 | Line 918 | namespace OpenMD {
918            // calculate the potential
919            pot_term = (ct_ij * c2ri - ctidotj * c3);
920            vterm = pref * pot_term;
921 <          *(idat.vpair) += vterm;
921 >          vpair += vterm;
922            epot +=  *(idat.sw)  * vterm;
923  
924            // calculate derivatives for the forces and torques
# Line 866 | Line 942 | namespace OpenMD {
942  
943          if (screeningMethod_ == DAMPED) {
944            // assemble the damping variables
945 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
946 <          erfcVal = res.first;
947 <          derfcVal = res.second;
945 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
946 >          //erfcVal = res.first;
947 >          //derfcVal = res.second;
948 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
949 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
950            c1 = erfcVal * riji;
951            c2 = (-derfcVal + c1) * riji;
952            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 894 | Line 972 | namespace OpenMD {
972                       qzz_i * (cz2 * c3 - c2ri) );
973          
974          vterm = pref * pot_term;
975 <        *(idat.vpair) += vterm;
975 >        vpair += vterm;
976          epot +=  *(idat.sw)  * vterm;
977  
978          // calculate the derivatives for the forces and torques
# Line 909 | Line 987 | namespace OpenMD {
987        }
988      }
989  
912 –    idat.pot[ELECTROSTATIC_FAMILY] += epot;
913 –    *(idat.f1) += dVdr;
990  
991 <    if (i_is_Dipole || i_is_Quadrupole)
992 <      *(idat.t1) -= cross(uz_i, duduz_i);
993 <    if (i_is_Quadrupole) {
994 <      *(idat.t1) -= cross(ux_i, dudux_i);
995 <      *(idat.t1) -= cross(uy_i, duduy_i);
996 <    }
997 <    
998 <    if (j_is_Dipole || j_is_Quadrupole)
999 <      *(idat.t2) -= cross(uz_j, duduz_j);
1000 <    if (j_is_Quadrupole) {
1001 <      *(idat.t2) -= cross(uz_j, dudux_j);
1002 <      *(idat.t2) -= cross(uz_j, duduy_j);
1003 <    }
991 >    if (!idat.excluded) {
992 >      *(idat.vpair) += vpair;
993 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
994 >      *(idat.f1) += dVdr;
995 >      
996 >      if (i_is_Dipole || i_is_Quadrupole)
997 >        *(idat.t1) -= cross(uz_i, duduz_i);
998 >      if (i_is_Quadrupole) {
999 >        *(idat.t1) -= cross(ux_i, dudux_i);
1000 >        *(idat.t1) -= cross(uy_i, duduy_i);
1001 >      }
1002 >      
1003 >      if (j_is_Dipole || j_is_Quadrupole)
1004 >        *(idat.t2) -= cross(uz_j, duduz_j);
1005 >      if (j_is_Quadrupole) {
1006 >        *(idat.t2) -= cross(uz_j, dudux_j);
1007 >        *(idat.t2) -= cross(uz_j, duduy_j);
1008 >      }
1009  
1010 <    return;
930 <  }  
1010 >    } else {
1011  
1012 <  void Electrostatic::calcSkipCorrection(InteractionData &idat) {
1012 >      // only accumulate the forces and torques resulting from the
1013 >      // indirect reaction field terms.
1014  
1015 <    if (!initialized_) initialize();
1016 <    
1017 <    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes->first];
937 <    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes->second];
938 <    
939 <    // logicals
940 <
941 <    bool i_is_Charge = data1.is_Charge;
942 <    bool i_is_Dipole = data1.is_Dipole;
943 <
944 <    bool j_is_Charge = data2.is_Charge;
945 <    bool j_is_Dipole = data2.is_Dipole;
946 <
947 <    RealType q_i, q_j;
948 <    
949 <    // The skippedCharge computation is needed by the real-space cutoff methods
950 <    // (i.e. shifted force and shifted potential)
951 <
952 <    if (i_is_Charge) {
953 <      q_i = data1.charge;
954 <      *(idat.skippedCharge2) += q_i;
955 <    }
956 <
957 <    if (j_is_Charge) {
958 <      q_j = data2.charge;
959 <      *(idat.skippedCharge1) += q_j;
960 <    }
961 <
962 <    // the rest of this function should only be necessary for reaction field.
963 <
964 <    if (summationMethod_ == esm_REACTION_FIELD) {
965 <      RealType riji, ri2, ri3;
966 <      RealType mu_i, ct_i;
967 <      RealType mu_j, ct_j;
968 <      RealType preVal, rfVal, vterm, dudr, pref, myPot(0.0);
969 <      Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
970 <
971 <      // some variables we'll need independent of electrostatic type:
1015 >      *(idat.vpair) += indirect_vpair;
1016 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
1017 >      *(idat.f1) += indirect_dVdr;
1018        
973 –      riji = 1.0 /  *(idat.rij) ;
974 –      rhat =  *(idat.d)  * riji;
975 –
976 –      if (i_is_Dipole) {
977 –        mu_i = data1.dipole_moment;
978 –        uz_i = idat.eFrame1->getColumn(2);      
979 –        ct_i = dot(uz_i, rhat);
980 –        duduz_i = V3Zero;
981 –      }
982 –            
983 –      if (j_is_Dipole) {
984 –        mu_j = data2.dipole_moment;
985 –        uz_j = idat.eFrame2->getColumn(2);      
986 –        ct_j = dot(uz_j, rhat);
987 –        duduz_j = V3Zero;
988 –      }
989 –    
990 –      if (i_is_Charge) {
991 –        if (j_is_Charge) {
992 –          preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
993 –          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij) ;
994 –          vterm = preVal * rfVal;
995 –          myPot +=  *(idat.sw)  * vterm;        
996 –          dudr  =  *(idat.sw)  * preVal * 2.0 * rfVal * riji;        
997 –          dVdr += dudr * rhat;
998 –        }
999 –        
1000 –        if (j_is_Dipole) {
1001 –          ri2 = riji * riji;
1002 –          ri3 = ri2 * riji;        
1003 –          pref =  *(idat.electroMult) * pre12_ * q_i * mu_j;
1004 –          vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
1005 –          myPot +=  *(idat.sw)  * vterm;        
1006 –          dVdr += - *(idat.sw)  * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1007 –          duduz_j += - *(idat.sw)  * pref * rhat * (ri2 - preRF2_ *  *(idat.rij) );
1008 –        }
1009 –      }
1010 –      if (i_is_Dipole) {
1011 –        if (j_is_Charge) {
1012 –          ri2 = riji * riji;
1013 –          ri3 = ri2 * riji;        
1014 –          pref =  *(idat.electroMult) * pre12_ * q_j * mu_i;
1015 –          vterm = - pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
1016 –          myPot +=  *(idat.sw)  * vterm;        
1017 –          dVdr +=  *(idat.sw)  * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);      
1018 –          duduz_i +=  *(idat.sw)  * pref * rhat * (ri2 - preRF2_ *  *(idat.rij));
1019 –        }
1020 –      }
1021 –      
1022 –      // accumulate the forces and torques resulting from the self term
1023 –      idat.pot[ELECTROSTATIC_FAMILY] += myPot;
1024 –      *(idat.f1) += dVdr;
1025 –      
1019        if (i_is_Dipole)
1020 <        *(idat.t1) -= cross(uz_i, duduz_i);
1020 >        *(idat.t1) -= cross(uz_i, indirect_duduz_i);
1021        if (j_is_Dipole)
1022 <        *(idat.t2) -= cross(uz_j, duduz_j);
1022 >        *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1023      }
1024 <  }
1024 >
1025 >
1026 >    return;
1027 >  }  
1028      
1029    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1030      RealType mu1, preVal, chg1, self;
1031      
1032      if (!initialized_) initialize();
1033 <    
1033 >
1034      ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1035    
1036      // logicals
1041 –
1037      bool i_is_Charge = data.is_Charge;
1038      bool i_is_Dipole = data.is_Dipole;
1039  
# Line 1046 | Line 1041 | namespace OpenMD {
1041        if (i_is_Dipole) {
1042          mu1 = data.dipole_moment;          
1043          preVal = pre22_ * preRF2_ * mu1 * mu1;
1044 <        sdat.pot[2] -= 0.5 * preVal;
1044 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal;
1045          
1046          // The self-correction term adds into the reaction field vector
1047          Vector3d uz_i = sdat.eFrame->getColumn(2);
# Line 1063 | Line 1058 | namespace OpenMD {
1058          } else {        
1059            self = - 0.5 * rcuti_ * chg1 * (chg1 +  *(sdat.skippedCharge)) * pre11_;
1060          }
1061 <        sdat.pot[ELECTROSTATIC_FAMILY] += self;
1061 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1062        }
1063      }
1064    }

Diff Legend

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