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

Comparing trunk/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1893 by gezelter, Wed Jun 19 17:19:07 2013 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 40 | Line 40
40   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43 + #ifdef IS_MPI
44 + #include <mpi.h>
45 + #endif
46 +
47   #include <stdio.h>
48   #include <string.h>
49  
50   #include <cmath>
51 + #include <numeric>
52   #include "nonbonded/Electrostatic.hpp"
53   #include "utils/simError.h"
54   #include "types/NonBondedInteractionType.hpp"
# Line 55 | Line 60
60   #include "utils/PhysicalConstants.hpp"
61   #include "math/erfc.hpp"
62   #include "math/SquareMatrix.hpp"
63 + #include "primitives/Molecule.hpp"
64 + #include "flucq/FluctuatingChargeForces.hpp"
65  
66   namespace OpenMD {
67    
68    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
62 –                                  forceField_(NULL), info_(NULL),
69                                    haveCutoffRadius_(false),
70                                    haveDampingAlpha_(false),
71                                    haveDielectric_(false),
72 <                                  haveElectroSplines_(false)
73 <  {}
72 >                                  haveElectroSplines_(false),
73 >                                  info_(NULL), forceField_(NULL)
74 >                                  
75 >  {
76 >    flucQ_ = new FluctuatingChargeForces(info_);
77 >  }
78    
79 +  void Electrostatic::setForceField(ForceField *ff) {
80 +    forceField_ = ff;
81 +    flucQ_->setForceField(forceField_);
82 +  }
83 +
84 +  void Electrostatic::setSimulatedAtomTypes(set<AtomType*> &simtypes) {
85 +    simTypes_ = simtypes;
86 +    flucQ_->setSimulatedAtomTypes(simTypes_);
87 +  }
88 +
89    void Electrostatic::initialize() {
90      
91      Globals* simParams_ = info_->getSimParams();
# Line 191 | Line 211 | namespace OpenMD {
211        simError();
212      }
213            
214 <    if (screeningMethod_ == DAMPED) {      
214 >    if (screeningMethod_ == DAMPED || summationMethod_ == esm_EWALD_FULL) {
215        if (!simParams_->haveDampingAlpha()) {
216          // first set a cutoff dependent alpha value
217          // we assume alpha depends linearly with rcut from 0 to 20.5 ang
218          dampingAlpha_ = 0.425 - cutoffRadius_* 0.02;
219 <        if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;
200 <        
219 >        if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;        
220          // throw warning
221          sprintf( painCave.errMsg,
222                   "Electrostatic::initialize: dampingAlpha was not specified in the\n"
# Line 213 | Line 232 | namespace OpenMD {
232        haveDampingAlpha_ = true;
233      }
234  
235 <    // find all of the Electrostatic atom Types:
236 <    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
237 <    ForceField::AtomTypeContainer::MapTypeIterator i;
238 <    AtomType* at;
235 >
236 >    Etypes.clear();
237 >    Etids.clear();
238 >    FQtypes.clear();
239 >    FQtids.clear();
240 >    ElectrostaticMap.clear();
241 >    Jij.clear();
242 >    nElectro_ = 0;
243 >    nFlucq_ = 0;
244 >
245 >    Etids.resize( forceField_->getNAtomType(), -1);
246 >    FQtids.resize( forceField_->getNAtomType(), -1);
247 >
248 >    set<AtomType*>::iterator at;
249 >    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {    
250 >      if ((*at)->isElectrostatic()) nElectro_++;
251 >      if ((*at)->isFluctuatingCharge()) nFlucq_++;
252 >    }
253      
254 <    for (at = atomTypes->beginType(i); at != NULL;
255 <         at = atomTypes->nextType(i)) {
256 <      
257 <      if (at->isElectrostatic())
225 <        addType(at);
254 >    Jij.resize(nFlucq_);
255 >
256 >    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
257 >      if ((*at)->isElectrostatic()) addType(*at);
258      }  
259      
260      if (summationMethod_ == esm_REACTION_FIELD) {
# Line 232 | Line 264 | namespace OpenMD {
264      
265      RealType b0c, b1c, b2c, b3c, b4c, b5c;
266      RealType db0c_1, db0c_2, db0c_3, db0c_4, db0c_5;
267 <    RealType a2, expTerm, invArootPi;
267 >    RealType a2, expTerm, invArootPi(0.0);
268      
269      RealType r = cutoffRadius_;
270      RealType r2 = r * r;
# Line 250 | Line 282 | namespace OpenMD {
282        b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2;
283        b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2;
284        b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2;
285 <      selfMult_ = b0c + a2 * invArootPi;
285 >      // Half the Smith self piece:
286 >      selfMult1_ = - a2 * invArootPi;
287 >      selfMult2_ = - 2.0 * a2 * a2 * invArootPi / 3.0;
288 >      selfMult4_ = - 4.0 * a2 * a2 * a2 * invArootPi / 5.0;
289      } else {
290        a2 = 0.0;
291        b0c = 1.0 / r;
# Line 259 | Line 294 | namespace OpenMD {
294        b3c = (5.0 * b2c) / r2;
295        b4c = (7.0 * b3c) / r2;
296        b5c = (9.0 * b4c) / r2;
297 <      selfMult_ = b0c;
297 >      selfMult1_ = 0.0;
298 >      selfMult2_ = 0.0;
299 >      selfMult4_ = 0.0;
300      }
301  
302      // higher derivatives of B_0 at the cutoff radius:
# Line 267 | Line 304 | namespace OpenMD {
304      db0c_2 =     -b1c + r2 * b2c;
305      db0c_3 =          3.0*r*b2c  - r2*r*b3c;
306      db0c_4 =          3.0*b2c  - 6.0*r2*b3c     + r2*r2*b4c;
307 <    db0c_5 =                    -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c;
271 <    
307 >    db0c_5 =                    -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c;  
308  
309 +    if (summationMethod_ != esm_EWALD_FULL) {
310 +      selfMult1_ -= b0c;
311 +      selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) /  3.0;
312 +      selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0;
313 +    }
314 +
315      // working variables for the splines:
316      RealType ri, ri2;
317      RealType b0, b1, b2, b3, b4, b5;
# Line 305 | Line 347 | namespace OpenMD {
347      vector<RealType> v31v, v32v;
348      vector<RealType> v41v, v42v, v43v;
349  
308 –    /*
309 –    vector<RealType> dv01v;
310 –    vector<RealType> dv11v;
311 –    vector<RealType> dv21v, dv22v;
312 –    vector<RealType> dv31v, dv32v;
313 –    vector<RealType> dv41v, dv42v, dv43v;
314 –    */
315 –
350      for (int i = 1; i < np_ + 1; i++) {
351        r = RealType(i) * dx;
352        rv.push_back(r);
# Line 398 | Line 432 | namespace OpenMD {
432          v11 = g - gc - rmRc*hc;
433          v21 = g*ri - gc*ric - rmRc*(hc - gc*ric)*ric;
434          v22 = h - g*ri - (hc - gc*ric) - rmRc*(sc - (hc - gc*ric)*ric);
435 <        v31 = (h-g*ri)*ri - (hc-g*ric)*ric - rmRc*(sc-2.0*(hc-gc*ric)*ric)*ric;
435 >        v31 = (h-g*ri)*ri - (hc-gc*ric)*ric - rmRc*(sc-2.0*(hc-gc*ric)*ric)*ric;
436          v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric)
437            - rmRc*(tc - 3.0*(sc-2.0*(hc-gc*ric)*ric)*ric);
438          v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2
# Line 455 | Line 489 | namespace OpenMD {
489          v11 = g - gc;
490          v21 = g*ri - gc*ric;
491          v22 = h - g*ri - (hc - gc*ric);
492 <        v31 = (h-g*ri)*ri - (hc-g*ric)*ric;
492 >        v31 = (h-g*ri)*ri - (hc-gc*ric)*ric;
493          v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
494          v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
495          v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;        
# Line 476 | Line 510 | namespace OpenMD {
510  
511        case esm_SWITCHING_FUNCTION:
512        case esm_HARD:
513 +      case esm_EWALD_FULL:
514  
515          v01 = f;
516          v11 = g;
# Line 515 | Line 550 | namespace OpenMD {
550          v11 = g - gc;
551          v21 = g*ri - gc*ric;
552          v22 = h - g*ri - (hc - gc*ric);
553 <        v31 = (h-g*ri)*ri - (hc-g*ric)*ric;
553 >        v31 = (h-g*ri)*ri - (hc-gc*ric)*ric;
554          v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
555          v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
556          v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;        
# Line 534 | Line 569 | namespace OpenMD {
569  
570          break;
571                  
537 –      case esm_EWALD_FULL:
572        case esm_EWALD_PME:
573        case esm_EWALD_SPME:
574        default :
# Line 563 | Line 597 | namespace OpenMD {
597        v41v.push_back(v41);
598        v42v.push_back(v42);
599        v43v.push_back(v43);
566 –      /*
567 –      dv01v.push_back(dv01);
568 –      dv11v.push_back(dv11);
569 –      dv21v.push_back(dv21);
570 –      dv22v.push_back(dv22);
571 –      dv31v.push_back(dv31);
572 –      dv32v.push_back(dv32);      
573 –      dv41v.push_back(dv41);
574 –      dv42v.push_back(dv42);
575 –      dv43v.push_back(dv43);
576 –      */
600      }
601  
602      // construct the spline structures and fill them with the values we've
# Line 598 | Line 621 | namespace OpenMD {
621      v43s = new CubicSpline();
622      v43s->addPoints(rv, v43v);
623  
601 –    /*
602 –    dv01s = new CubicSpline();
603 –    dv01s->addPoints(rv, dv01v);
604 –    dv11s = new CubicSpline();
605 –    dv11s->addPoints(rv, dv11v);
606 –    dv21s = new CubicSpline();
607 –    dv21s->addPoints(rv, dv21v);
608 –    dv22s = new CubicSpline();
609 –    dv22s->addPoints(rv, dv22v);
610 –    dv31s = new CubicSpline();
611 –    dv31s->addPoints(rv, dv31v);
612 –    dv32s = new CubicSpline();
613 –    dv32s->addPoints(rv, dv32v);
614 –    dv41s = new CubicSpline();
615 –    dv41s->addPoints(rv, dv41v);
616 –    dv42s = new CubicSpline();
617 –    dv42s->addPoints(rv, dv42v);
618 –    dv43s = new CubicSpline();
619 –    dv43s->addPoints(rv, dv43v);
620 –    */
621 –
624      haveElectroSplines_ = true;
625  
626      initialized_ = true;
627    }
628        
629    void Electrostatic::addType(AtomType* atomType){
630 <
630 >    
631      ElectrostaticAtomData electrostaticAtomData;
632      electrostaticAtomData.is_Charge = false;
633      electrostaticAtomData.is_Dipole = false;
# Line 661 | Line 663 | namespace OpenMD {
663        electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
664      }
665  
666 <    pair<map<int,AtomType*>::iterator,bool> ret;    
667 <    ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
668 <                                                        atomType) );
666 >    int atid = atomType->getIdent();
667 >    int etid = Etypes.size();
668 >    int fqtid = FQtypes.size();
669 >
670 >    pair<set<int>::iterator,bool> ret;    
671 >    ret = Etypes.insert( atid );
672      if (ret.second == false) {
673        sprintf( painCave.errMsg,
674                 "Electrostatic already had a previous entry with ident %d\n",
675 <               atomType->getIdent() );
675 >               atid);
676        painCave.severity = OPENMD_INFO;
677        painCave.isFatal = 0;
678        simError();        
679      }
680      
681 <    ElectrostaticMap[atomType] = electrostaticAtomData;  
681 >    Etids[ atid ] = etid;
682 >    ElectrostaticMap.push_back(electrostaticAtomData);
683  
684 <    // Now, iterate over all known types and add to the mixing map:
685 <    
686 <    map<AtomType*, ElectrostaticAtomData>::iterator it;
687 <    for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
688 <      AtomType* atype2 = (*it).first;
689 <      ElectrostaticAtomData eaData2 = (*it).second;
690 <      if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) {
691 <        
684 >    if (electrostaticAtomData.is_Fluctuating) {
685 >      ret = FQtypes.insert( atid );
686 >      if (ret.second == false) {
687 >        sprintf( painCave.errMsg,
688 >                 "Electrostatic already had a previous fluctuating charge entry with ident %d\n",
689 >                 atid );
690 >        painCave.severity = OPENMD_INFO;
691 >        painCave.isFatal = 0;
692 >        simError();        
693 >      }
694 >      FQtids[atid] = fqtid;
695 >      Jij[fqtid].resize(nFlucq_);
696 >
697 >      // Now, iterate over all known fluctuating and add to the
698 >      // coulomb integral map:
699 >      
700 >      std::set<int>::iterator it;
701 >      for( it = FQtypes.begin(); it != FQtypes.end(); ++it) {    
702 >        int etid2 = Etids[ (*it) ];
703 >        int fqtid2 = FQtids[ (*it) ];
704 >        ElectrostaticAtomData eaData2 = ElectrostaticMap[ etid2 ];
705          RealType a = electrostaticAtomData.slaterZeta;
706          RealType b = eaData2.slaterZeta;
707          int m = electrostaticAtomData.slaterN;
708          int n = eaData2.slaterN;
709 <
709 >        
710          // Create the spline of the coulombic integral for s-type
711          // Slater orbitals.  Add a 2 angstrom safety window to deal
712          // with cutoffGroups that have charged atoms longer than the
713          // cutoffRadius away from each other.
714 <
714 >        
715          RealType rval;
716          RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
717          vector<RealType> rvals;
# Line 709 | Line 728 | namespace OpenMD {
728          
729          CubicSpline* J = new CubicSpline();
730          J->addPoints(rvals, Jvals);
731 <        
732 <        pair<AtomType*, AtomType*> key1, key2;
733 <        key1 = make_pair(atomType, atype2);
734 <        key2 = make_pair(atype2, atomType);
735 <        
717 <        Jij[key1] = J;
718 <        Jij[key2] = J;
719 <      }
720 <    }
721 <
731 >        Jij[fqtid][fqtid2] = J;
732 >        Jij[fqtid2].resize( nFlucq_ );
733 >        Jij[fqtid2][fqtid] = J;
734 >      }      
735 >    }      
736      return;
737    }
738    
# Line 746 | Line 760 | namespace OpenMD {
760  
761      if (!initialized_) initialize();
762      
763 <    data1 = ElectrostaticMap[idat.atypes.first];
764 <    data2 = ElectrostaticMap[idat.atypes.second];
763 >    data1 = ElectrostaticMap[Etids[idat.atid1]];
764 >    data2 = ElectrostaticMap[Etids[idat.atid2]];
765  
766      U = 0.0;  // Potential
767      F.zero();  // Force
# Line 755 | Line 769 | namespace OpenMD {
769      Tb.zero(); // Torque on site b
770      Ea.zero(); // Electric field at site a
771      Eb.zero(); // Electric field at site b
772 +    Pa = 0.0;  // Site potential at site a
773 +    Pb = 0.0;  // Site potential at site b
774      dUdCa = 0.0; // fluctuating charge force at site a
775      dUdCb = 0.0; // fluctuating charge force at site a
776      
# Line 767 | Line 783 | namespace OpenMD {
783      // Excluded potential that is still computed for fluctuating charges
784      excluded_Pot= 0.0;
785  
770 –
786      // some variables we'll need independent of electrostatic type:
787  
788      ri = 1.0 /  *(idat.rij);
# Line 830 | Line 845 | namespace OpenMD {
845        if (idat.excluded) {
846          *(idat.skippedCharge2) += C_a;
847        } else {
848 <        // only do the field if we're not excluded:
848 >        // only do the field and site potentials if we're not excluded:
849          Eb -= C_a *  pre11_ * dv01 * rhat;
850 +        Pb += C_a *  pre11_ * v01;
851        }
852      }
853      
# Line 839 | Line 855 | namespace OpenMD {
855        D_a = *(idat.dipole1);
856        rdDa = dot(rhat, D_a);
857        rxDa = cross(rhat, D_a);
858 <      if (!idat.excluded)
858 >      if (!idat.excluded) {
859          Eb -=  pre12_ * ((dv11-v11or) * rdDa * rhat + v11or * D_a);
860 +        Pb +=  pre12_ * v11 * rdDa;
861 +      }
862 +
863      }
864      
865      if (a_is_Quadrupole) {
# Line 850 | Line 869 | namespace OpenMD {
869        rQa = rhat * Q_a;
870        rdQar = dot(rhat, Qar);
871        rxQar = cross(rhat, Qar);
872 <      if (!idat.excluded)
872 >      if (!idat.excluded) {
873          Eb -= pre14_ * (trQa * rhat * dv21 + 2.0 * Qar * v22or
874                          + rdQar * rhat * (dv22 - 2.0*v22or));
875 +        Pb += pre14_ * (v21 * trQa + v22 * rdQar);
876 +      }
877      }
878      
879      if (b_is_Charge) {
# Line 866 | Line 887 | namespace OpenMD {
887        } else {
888          // only do the field if we're not excluded:
889          Ea += C_b *  pre11_ * dv01 * rhat;
890 +        Pa += C_b *  pre11_ * v01;
891 +
892        }
893      }
894      
# Line 873 | Line 896 | namespace OpenMD {
896        D_b = *(idat.dipole2);
897        rdDb = dot(rhat, D_b);
898        rxDb = cross(rhat, D_b);
899 <      if (!idat.excluded)
899 >      if (!idat.excluded) {
900          Ea += pre12_ * ((dv11-v11or) * rdDb * rhat + v11or * D_b);
901 +        Pa += pre12_ * v11 * rdDb;
902 +      }
903      }
904      
905      if (b_is_Quadrupole) {
# Line 884 | Line 909 | namespace OpenMD {
909        rQb = rhat * Q_b;
910        rdQbr = dot(rhat, Qbr);
911        rxQbr = cross(rhat, Qbr);
912 <      if (!idat.excluded)
912 >      if (!idat.excluded) {
913          Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or
914                          + rdQbr * rhat * (dv22 - 2.0*v22or));
915 +        Pa += pre14_ * (v21 * trQb + v22 * rdQbr);
916 +      }
917      }
918 <    
918 >        
919 >
920      if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
921 <      J = Jij[idat.atypes];
921 >      J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
922      }    
923 <    
923 >
924      if (a_is_Charge) {    
925        
926        if (b_is_Charge) {
927          pref =  pre11_ * *(idat.electroMult);      
928          U  += C_a * C_b * pref * v01;
929          F  += C_a * C_b * pref * dv01 * rhat;
930 <        
930 >
931          // If this is an excluded pair, there are still indirect
932          // interactions via the reaction field we must worry about:
933  
# Line 908 | Line 936 | namespace OpenMD {
936            indirect_Pot += rfContrib;
937            indirect_F   += rfContrib * 2.0 * ri * rhat;
938          }
939 <        
939 >
940          // Fluctuating charge forces are handled via Coulomb integrals
941          // for excluded pairs (i.e. those connected via bonds) and
942          // with the standard charge-charge interaction otherwise.
943  
944 <        if (idat.excluded) {          
944 >        if (idat.excluded) {
945            if (a_is_Fluctuating || b_is_Fluctuating) {
946              coulInt = J->getValueAt( *(idat.rij) );
947 <            if (a_is_Fluctuating)  dUdCa += coulInt * C_b;
948 <            if (b_is_Fluctuating)  dUdCb += coulInt * C_a;
949 <            excluded_Pot += C_a * C_b * coulInt;
922 <          }          
947 >            if (a_is_Fluctuating) dUdCa += C_b * coulInt;
948 >            if (b_is_Fluctuating) dUdCb += C_a * coulInt;          
949 >          }
950          } else {
951            if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
952 <          if (a_is_Fluctuating) dUdCb += C_a * pref * v01;
953 <        }
952 >          if (b_is_Fluctuating) dUdCb += C_a * pref * v01;
953 >        }              
954        }
955  
956        if (b_is_Dipole) {
# Line 989 | Line 1016 | namespace OpenMD {
1016          F  -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat;
1017          Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
1018          Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
992 –
1019          // Even if we excluded this pair from direct interactions, we
1020          // still have the reaction-field-mediated dipole-dipole
1021          // interaction:
# Line 1049 | Line 1075 | namespace OpenMD {
1075          trQaQb = QaQb.trace();
1076          rQaQb = rhat * QaQb;
1077          QaQbr = QaQb * rhat;
1078 <        QaxQb = cross(Q_a, Q_b);
1078 >        QaxQb = mCross(Q_a, Q_b);
1079          rQaQbr = dot(rQa, Qbr);
1080          rQaxQbr = cross(rQa, Qbr);
1081          
# Line 1080 | Line 1106 | namespace OpenMD {
1106          //             + 4.0 * cross(rhat, QbQar)
1107  
1108          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1083 –
1084 –        //  cerr << " tsum = " << Ta + Tb - cross(  *(idat.d) , F ) << "\n";
1109        }
1110      }
1111  
# Line 1090 | Line 1114 | namespace OpenMD {
1114        *(idat.eField2) += Eb * *(idat.electroMult);
1115      }
1116  
1117 +    if (idat.doSitePotential) {
1118 +      *(idat.sPot1) += Pa * *(idat.electroMult);
1119 +      *(idat.sPot2) += Pb * *(idat.electroMult);
1120 +    }
1121 +
1122      if (a_is_Fluctuating) *(idat.dVdFQ1) += dUdCa * *(idat.sw);
1123      if (b_is_Fluctuating) *(idat.dVdFQ2) += dUdCb * *(idat.sw);
1124  
# Line 1128 | Line 1157 | namespace OpenMD {
1157  
1158      if (!initialized_) initialize();
1159  
1160 <    ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1160 >    ElectrostaticAtomData data = ElectrostaticMap[Etids[sdat.atid]];
1161      
1162      // logicals
1163      bool i_is_Charge = data.is_Charge;
1164      bool i_is_Dipole = data.is_Dipole;
1165 +    bool i_is_Quadrupole = data.is_Quadrupole;
1166      bool i_is_Fluctuating = data.is_Fluctuating;
1167      RealType C_a = data.fixedCharge;  
1168 <    RealType self, preVal, DadDa;
1169 <    
1168 >    RealType self(0.0), preVal, DdD(0.0), trQ, trQQ;
1169 >
1170 >    if (i_is_Dipole) {
1171 >      DdD = data.dipole.lengthSquare();
1172 >    }
1173 >        
1174      if (i_is_Fluctuating) {
1175        C_a += *(sdat.flucQ);
1176 <      // dVdFQ is really a force, so this is negative the derivative
1177 <      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1178 <      (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1179 <        (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1176 >
1177 >      flucQ_->getSelfInteraction(sdat.atid, *(sdat.flucQ),  
1178 >                                 (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY],
1179 >                                 *(sdat.flucQfrc) );
1180 >
1181      }
1182  
1183      switch (summationMethod_) {
# Line 1157 | Line 1192 | namespace OpenMD {
1192        }
1193  
1194        if (i_is_Dipole) {
1195 <        DadDa = data.dipole.lengthSquare();
1161 <        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DadDa;
1195 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DdD;
1196        }
1197        
1198        break;
1199        
1200      case esm_SHIFTED_FORCE:
1201      case esm_SHIFTED_POTENTIAL:
1202 <      if (i_is_Charge) {
1203 <        self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
1204 <        (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1202 >    case esm_TAYLOR_SHIFTED:
1203 >    case esm_EWALD_FULL:
1204 >      if (i_is_Charge)
1205 >        self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge));      
1206 >      if (i_is_Dipole)
1207 >        self += selfMult2_ * pre22_ * DdD;      
1208 >      if (i_is_Quadrupole) {
1209 >        trQ = data.quadrupole.trace();
1210 >        trQQ = (data.quadrupole * data.quadrupole).trace();
1211 >        self += selfMult4_ * pre44_ * (2.0*trQQ + trQ*trQ);
1212 >        if (i_is_Charge)
1213 >          self -= selfMult2_ * pre14_ * 2.0 * C_a * trQ;
1214        }
1215 +      (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;      
1216        break;
1217      default:
1218        break;
# Line 1182 | Line 1226 | namespace OpenMD {
1226      // cases.
1227      return 12.0;
1228    }
1229 +
1230 +
1231 +  void Electrostatic::ReciprocalSpaceSum(RealType& pot) {
1232 +    
1233 +    RealType kPot = 0.0;
1234 +    RealType kVir = 0.0;
1235 +    
1236 +    const RealType mPoleConverter = 0.20819434; // converts from the
1237 +                                                // internal units of
1238 +                                                // Debye (for dipoles)
1239 +                                                // or Debye-angstroms
1240 +                                                // (for quadrupoles) to
1241 +                                                // electron angstroms or
1242 +                                                // electron-angstroms^2
1243 +    
1244 +    const RealType eConverter = 332.0637778; // convert the
1245 +                                             // Charge-Charge
1246 +                                             // electrostatic
1247 +                                             // interactions into kcal /
1248 +                                             // mol assuming distances
1249 +                                             // are measured in
1250 +                                             // angstroms.
1251 +
1252 +    Mat3x3d hmat = info_->getSnapshotManager()->getCurrentSnapshot()->getHmat();
1253 +    Vector3d box = hmat.diagonals();
1254 +    RealType boxMax = box.max();
1255 +    
1256 +    //int kMax = int(2.0 * M_PI / (pow(dampingAlpha_,2)*cutoffRadius_ * boxMax) );
1257 +    int kMax = 7;
1258 +    int kSqMax = kMax*kMax + 2;
1259 +    
1260 +    int kLimit = kMax+1;
1261 +    int kLim2 = 2*kMax+1;
1262 +    int kSqLim = kSqMax;
1263 +    
1264 +    vector<RealType> AK(kSqLim+1, 0.0);
1265 +    RealType xcl = 2.0 * M_PI / box.x();
1266 +    RealType ycl = 2.0 * M_PI / box.y();
1267 +    RealType zcl = 2.0 * M_PI / box.z();
1268 +    RealType rcl = 2.0 * M_PI / boxMax;
1269 +    RealType rvol = 2.0 * M_PI /(box.x() * box.y() * box.z());
1270 +    
1271 +    if(dampingAlpha_ < 1.0e-12) return;
1272 +    
1273 +    RealType ralph = -0.25/pow(dampingAlpha_,2);
1274 +    
1275 +    // Calculate and store exponential factors  
1276 +    
1277 +    vector<vector<RealType> > elc;
1278 +    vector<vector<RealType> > emc;
1279 +    vector<vector<RealType> > enc;
1280 +    vector<vector<RealType> > els;
1281 +    vector<vector<RealType> > ems;
1282 +    vector<vector<RealType> > ens;
1283 +    
1284 +    int nMax = info_->getNAtoms();
1285 +    
1286 +    elc.resize(kLimit+1);
1287 +    emc.resize(kLimit+1);
1288 +    enc.resize(kLimit+1);
1289 +    els.resize(kLimit+1);
1290 +    ems.resize(kLimit+1);
1291 +    ens.resize(kLimit+1);
1292 +
1293 +    for (int j = 0; j < kLimit+1; j++) {
1294 +      elc[j].resize(nMax);
1295 +      emc[j].resize(nMax);
1296 +      enc[j].resize(nMax);
1297 +      els[j].resize(nMax);
1298 +      ems[j].resize(nMax);
1299 +      ens[j].resize(nMax);
1300 +    }
1301 +    
1302 +    Vector3d t( 2.0 * M_PI );
1303 +    t.Vdiv(t, box);
1304 +
1305 +    SimInfo::MoleculeIterator mi;
1306 +    Molecule::AtomIterator ai;
1307 +    int i;
1308 +    Vector3d r;
1309 +    Vector3d tt;
1310 +    
1311 +    for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1312 +         mol = info_->nextMolecule(mi)) {
1313 +      for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1314 +          atom = mol->nextAtom(ai)) {  
1315 +        
1316 +        i = atom->getLocalIndex();
1317 +        r = atom->getPos();
1318 +        info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r);
1319 +        
1320 +        tt.Vmul(t, r);
1321 +
1322 +        elc[1][i] = 1.0;
1323 +        emc[1][i] = 1.0;
1324 +        enc[1][i] = 1.0;
1325 +        els[1][i] = 0.0;
1326 +        ems[1][i] = 0.0;
1327 +        ens[1][i] = 0.0;
1328 +
1329 +        elc[2][i] = cos(tt.x());
1330 +        emc[2][i] = cos(tt.y());
1331 +        enc[2][i] = cos(tt.z());
1332 +        els[2][i] = sin(tt.x());
1333 +        ems[2][i] = sin(tt.y());
1334 +        ens[2][i] = sin(tt.z());
1335 +        
1336 +        for(int l = 3; l <= kLimit; l++) {
1337 +          elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i];
1338 +          emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i];
1339 +          enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i];
1340 +          els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i];
1341 +          ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i];
1342 +          ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i];
1343 +        }
1344 +      }
1345 +    }
1346 +    
1347 +    // Calculate and store AK coefficients:
1348 +    
1349 +    RealType eksq = 1.0;
1350 +    RealType expf = 0.0;
1351 +    if (ralph < 0.0) expf = exp(ralph*rcl*rcl);
1352 +    for (i = 1; i <= kSqLim; i++) {
1353 +      RealType rksq = float(i)*rcl*rcl;
1354 +      eksq = expf*eksq;
1355 +      AK[i] = eConverter * eksq/rksq;
1356 +    }
1357 +    
1358 +    /*
1359 +     * Loop over all k vectors k = 2 pi (ll/Lx, mm/Ly, nn/Lz)
1360 +     * the values of ll, mm and nn are selected so that the symmetry of
1361 +     * reciprocal lattice is taken into account i.e. the following
1362 +     * rules apply.
1363 +     *
1364 +     * ll ranges over the values 0 to kMax only.
1365 +     *
1366 +     * mm ranges over 0 to kMax when ll=0 and over
1367 +     *            -kMax to kMax otherwise.
1368 +     * nn ranges over 1 to kMax when ll=mm=0 and over
1369 +     *            -kMax to kMax otherwise.
1370 +     *
1371 +     * Hence the result of the summation must be doubled at the end.    
1372 +     */
1373 +    
1374 +    std::vector<RealType> clm(nMax, 0.0);
1375 +    std::vector<RealType> slm(nMax, 0.0);
1376 +    std::vector<RealType> ckr(nMax, 0.0);
1377 +    std::vector<RealType> skr(nMax, 0.0);
1378 +    std::vector<RealType> ckc(nMax, 0.0);
1379 +    std::vector<RealType> cks(nMax, 0.0);
1380 +    std::vector<RealType> dkc(nMax, 0.0);
1381 +    std::vector<RealType> dks(nMax, 0.0);
1382 +    std::vector<RealType> qkc(nMax, 0.0);
1383 +    std::vector<RealType> qks(nMax, 0.0);
1384 +    std::vector<Vector3d> dxk(nMax, V3Zero);
1385 +    std::vector<Vector3d> qxk(nMax, V3Zero);
1386 +    RealType rl, rm, rn;
1387 +    Vector3d kVec;
1388 +    Vector3d Qk;
1389 +    Mat3x3d k2;
1390 +    RealType ckcs, ckss, dkcs, dkss, qkcs, qkss;
1391 +    int atid;
1392 +    ElectrostaticAtomData data;
1393 +    RealType C, dk, qk;
1394 +    Vector3d D;
1395 +    Mat3x3d  Q;
1396 +
1397 +    int mMin = kLimit;
1398 +    int nMin = kLimit + 1;
1399 +    for (int l = 1; l <= kLimit; l++) {
1400 +      int ll = l - 1;
1401 +      rl = xcl * float(ll);
1402 +      for (int mmm = mMin; mmm <= kLim2; mmm++) {
1403 +        int mm = mmm - kLimit;
1404 +        int m = abs(mm) + 1;
1405 +        rm = ycl * float(mm);
1406 +        // Set temporary products of exponential terms
1407 +        for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1408 +             mol = info_->nextMolecule(mi)) {
1409 +          for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1410 +              atom = mol->nextAtom(ai)) {
1411 +            
1412 +            i = atom->getLocalIndex();
1413 +            if(mm < 0) {
1414 +              clm[i]=elc[l][i]*emc[m][i]+els[l][i]*ems[m][i];
1415 +              slm[i]=els[l][i]*emc[m][i]-ems[m][i]*elc[l][i];
1416 +            } else {
1417 +              clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i];
1418 +              slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i];
1419 +            }
1420 +          }
1421 +        }
1422 +        for (int nnn = nMin; nnn <= kLim2; nnn++) {
1423 +          int nn = nnn - kLimit;          
1424 +          int n = abs(nn) + 1;
1425 +          rn = zcl * float(nn);
1426 +          // Test on magnitude of k vector:
1427 +          int kk=ll*ll + mm*mm + nn*nn;
1428 +          if(kk <= kSqLim) {
1429 +            kVec = Vector3d(rl, rm, rn);
1430 +            k2 = outProduct(kVec, kVec);
1431 +            // Calculate exp(ikr) terms
1432 +            for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1433 +                 mol = info_->nextMolecule(mi)) {
1434 +              for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1435 +                  atom = mol->nextAtom(ai)) {
1436 +                i = atom->getLocalIndex();
1437 +                
1438 +                if (nn < 0) {
1439 +                  ckr[i]=clm[i]*enc[n][i]+slm[i]*ens[n][i];
1440 +                  skr[i]=slm[i]*enc[n][i]-clm[i]*ens[n][i];
1441 +
1442 +                } else {
1443 +                  ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i];
1444 +                  skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i];
1445 +                }
1446 +              }
1447 +            }
1448 +            
1449 +            // Calculate scalar and vector products for each site:
1450 +            
1451 +            for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1452 +                 mol = info_->nextMolecule(mi)) {
1453 +              for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1454 +                  atom = mol->nextAtom(ai)) {
1455 +                i = atom->getLocalIndex();
1456 +                int atid = atom->getAtomType()->getIdent();
1457 +                data = ElectrostaticMap[Etids[atid]];
1458 +                              
1459 +                if (data.is_Charge) {
1460 +                  C = data.fixedCharge;
1461 +                  if (atom->isFluctuatingCharge()) C += atom->getFlucQPos();
1462 +                  ckc[i] = C * ckr[i];
1463 +                  cks[i] = C * skr[i];
1464 +                }
1465 +                
1466 +                if (data.is_Dipole) {
1467 +                  D = atom->getDipole() * mPoleConverter;
1468 +                  dk = dot(D, kVec);
1469 +                  dxk[i] = cross(D, kVec);
1470 +                  dkc[i] = dk * ckr[i];
1471 +                  dks[i] = dk * skr[i];
1472 +                }
1473 +                if (data.is_Quadrupole) {
1474 +                  Q = atom->getQuadrupole() * mPoleConverter;
1475 +                  Qk = Q * kVec;                  
1476 +                  qk = dot(kVec, Qk);
1477 +                  qxk[i] = -cross(kVec, Qk);
1478 +                  qkc[i] = qk * ckr[i];
1479 +                  qks[i] = qk * skr[i];
1480 +                }              
1481 +              }
1482 +            }
1483 +
1484 +            // calculate vector sums
1485 +            
1486 +            ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1487 +            ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1488 +            dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1489 +            dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1490 +            qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1491 +            qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1492 +            
1493 + #ifdef IS_MPI
1494 +            MPI_Allreduce(MPI_IN_PLACE, &ckcs, 1, MPI_REALTYPE,
1495 +                          MPI_SUM, MPI_COMM_WORLD);
1496 +            MPI_Allreduce(MPI_IN_PLACE, &ckss, 1, MPI_REALTYPE,
1497 +                          MPI_SUM, MPI_COMM_WORLD);
1498 +            MPI_Allreduce(MPI_IN_PLACE, &dkcs, 1, MPI_REALTYPE,
1499 +                          MPI_SUM, MPI_COMM_WORLD);
1500 +            MPI_Allreduce(MPI_IN_PLACE, &dkss, 1, MPI_REALTYPE,
1501 +                          MPI_SUM, MPI_COMM_WORLD);
1502 +            MPI_Allreduce(MPI_IN_PLACE, &qkcs, 1, MPI_REALTYPE,
1503 +                          MPI_SUM, MPI_COMM_WORLD);
1504 +            MPI_Allreduce(MPI_IN_PLACE, &qkss, 1, MPI_REALTYPE,
1505 +                          MPI_SUM, MPI_COMM_WORLD);
1506 + #endif        
1507 +            
1508 +            // Accumulate potential energy and virial contribution:
1509 +
1510 +            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1511 +                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs));
1512 +
1513 +            kVir += 2.0 * rvol  * AK[kk]*(ckcs*ckcs+ckss*ckss
1514 +                                          +4.0*(ckss*dkcs-ckcs*dkss)
1515 +                                          +3.0*(dkcs*dkcs+dkss*dkss)
1516 +                                          -6.0*(ckss*qkss+ckcs*qkcs)
1517 +                                          +8.0*(dkss*qkcs-dkcs*qkss)
1518 +                                          +5.0*(qkss*qkss+qkcs*qkcs));
1519 +            
1520 +            // Calculate force and torque for each site:
1521 +            
1522 +            for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1523 +                 mol = info_->nextMolecule(mi)) {
1524 +              for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1525 +                  atom = mol->nextAtom(ai)) {
1526 +                
1527 +                i = atom->getLocalIndex();
1528 +                atid = atom->getAtomType()->getIdent();
1529 +                data = ElectrostaticMap[Etids[atid]];
1530 +
1531 +                RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs)
1532 +                                     - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss));
1533 +                RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
1534 +                                         -ckr[i]*(ckss+dkcs-qkss));
1535 +                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
1536 +                                            +skr[i]*(ckss+dkcs-qkss));
1537 +              
1538 +                atom->addFrc( 4.0 * rvol * qfrc * kVec );
1539 +
1540 +                if (atom->isFluctuatingCharge()) {
1541 +                  atom->addFlucQFrc( - 2.0 * rvol * qtrq2 );
1542 +                }
1543 +                  
1544 +                if (data.is_Dipole) {
1545 +                  atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] );
1546 +                }
1547 +                if (data.is_Quadrupole) {
1548 +                  atom->addTrq( 4.0 * rvol * qtrq2 * qxk[i] );
1549 +                }
1550 +              }
1551 +            }
1552 +          }
1553 +        }
1554 +        nMin = 1;
1555 +      }
1556 +      mMin = 1;
1557 +    }
1558 +    pot += kPot;  
1559 +  }
1560 +
1561 +  void Electrostatic::getSitePotentials(Atom* a1, Atom* a2, bool excluded,
1562 +                                        RealType &spot1, RealType &spot2) {
1563 +
1564 +    if (!initialized_) {
1565 +      cerr << "initializing\n";
1566 +      initialize();
1567 +      cerr << "done\n";
1568 +    }
1569 +
1570 +    const RealType mPoleConverter = 0.20819434;
1571 +
1572 +    AtomType* atype1 = a1->getAtomType();
1573 +    AtomType* atype2 = a2->getAtomType();
1574 +    int atid1 = atype1->getIdent();
1575 +    int atid2 = atype2->getIdent();
1576 +    data1 = ElectrostaticMap[Etids[atid1]];
1577 +    data2 = ElectrostaticMap[Etids[atid2]];
1578 +
1579 +    Pa = 0.0;  // Site potential at site a
1580 +    Pb = 0.0;  // Site potential at site b
1581 +
1582 +    Vector3d d = a2->getPos() - a1->getPos();
1583 +    info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(d);
1584 +    RealType rij = d.length();
1585 +    // some variables we'll need independent of electrostatic type:
1586 +
1587 +    RealType ri = 1.0 /  rij;
1588 +    rhat = d  * ri;
1589 +      
1590 +
1591 +    if ((rij >= cutoffRadius_) || excluded) {
1592 +      spot1 = 0.0;
1593 +      spot2 = 0.0;
1594 +      return;
1595 +    }
1596 +
1597 +    // logicals
1598 +
1599 +    a_is_Charge = data1.is_Charge;
1600 +    a_is_Dipole = data1.is_Dipole;
1601 +    a_is_Quadrupole = data1.is_Quadrupole;
1602 +    a_is_Fluctuating = data1.is_Fluctuating;
1603 +
1604 +    b_is_Charge = data2.is_Charge;
1605 +    b_is_Dipole = data2.is_Dipole;
1606 +    b_is_Quadrupole = data2.is_Quadrupole;
1607 +    b_is_Fluctuating = data2.is_Fluctuating;
1608 +
1609 +    // Obtain all of the required radial function values from the
1610 +    // spline structures:
1611 +    
1612 +
1613 +    if (a_is_Charge || b_is_Charge) {
1614 +      v01 = v01s->getValueAt(rij);
1615 +    }
1616 +    if (a_is_Dipole || b_is_Dipole) {
1617 +      v11 = v11s->getValueAt(rij);
1618 +      v11or = ri * v11;
1619 +    }
1620 +    if (a_is_Quadrupole || b_is_Quadrupole) {
1621 +      v21 = v21s->getValueAt(rij);
1622 +      v22 = v22s->getValueAt(rij);
1623 +      v22or = ri * v22;
1624 +    }      
1625 +
1626 +    if (a_is_Charge) {
1627 +      C_a = data1.fixedCharge;
1628 +      
1629 +      if (a_is_Fluctuating) {
1630 +        C_a += a1->getFlucQPos();
1631 +      }
1632 +      
1633 +      Pb += C_a *  pre11_ * v01;      
1634 +    }
1635 +    
1636 +    if (a_is_Dipole) {
1637 +      D_a = a1->getDipole() * mPoleConverter;
1638 +      rdDa = dot(rhat, D_a);
1639 +      Pb +=  pre12_ * v11 * rdDa;      
1640 +    }
1641 +    
1642 +    if (a_is_Quadrupole) {
1643 +      Q_a = a1->getQuadrupole() * mPoleConverter;
1644 +      trQa =  Q_a.trace();
1645 +      Qar =   Q_a * rhat;
1646 +      rdQar = dot(rhat, Qar);
1647 +      Pb += pre14_ * (v21 * trQa + v22 * rdQar);      
1648 +    }
1649 +    
1650 +    if (b_is_Charge) {
1651 +      C_b = data2.fixedCharge;
1652 +      
1653 +      if (b_is_Fluctuating)
1654 +        C_b += a2->getFlucQPos();
1655 +      
1656 +      Pa += C_b *  pre11_ * v01;
1657 +    }
1658 +    
1659 +    if (b_is_Dipole) {
1660 +      D_b = a2->getDipole() * mPoleConverter;
1661 +      rdDb = dot(rhat, D_b);
1662 +      Pa += pre12_ * v11 * rdDb;
1663 +    }
1664 +    
1665 +    if (b_is_Quadrupole) {
1666 +      Q_a = a2->getQuadrupole() * mPoleConverter;
1667 +      trQb =  Q_b.trace();
1668 +      Qbr =   Q_b * rhat;
1669 +      rdQbr = dot(rhat, Qbr);
1670 +      Pa += pre14_ * (v21 * trQb + v22 * rdQbr);      
1671 +    }
1672 +
1673 +    spot1 = Pa;
1674 +    spot2 = Pb;
1675 +  }
1676   }

Diff Legend

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