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 1528 by gezelter, Fri Dec 17 20:11:05 2010 UTC vs.
Revision 1601 by gezelter, Thu Aug 4 20:04:35 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 47 | Line 47
47   #include "utils/simError.h"
48   #include "types/NonBondedInteractionType.hpp"
49   #include "types/DirectionalAtomType.hpp"
50 + #include "io/Globals.hpp"
51  
51 –
52   namespace OpenMD {
53    
54    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
55 <                                  forceField_(NULL), info_(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  
66      summationMap_["HARD"]               = esm_HARD;
67      summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
# Line 95 | Line 102 | namespace OpenMD {
102      screeningMethod_ = UNDAMPED;
103      dielectric_ = 1.0;
104      one_third_ = 1.0 / 3.0;
98 –    haveCutoffRadius_ = false;
99 –    haveDampingAlpha_ = false;
100 –    haveDielectric_ = false;  
101 –    haveElectroSpline_ = false;
105    
106      // check the summation method:
107      if (simParams_->haveElectrostaticSummationMethod()) {
# Line 111 | 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 405 | Line 408 | namespace OpenMD {
408      return;
409    }
410    
411 <  void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
412 <                                                    RealType theRSW ) {
410 <    cutoffRadius_ = theECR;
411 >  void Electrostatic::setCutoffRadius( RealType rCut ) {
412 >    cutoffRadius_ = rCut;
413      rrf_ = cutoffRadius_;
412 –    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 427 | Line 432 | namespace OpenMD {
432      haveDielectric_ = true;
433    }
434  
435 <  void Electrostatic::calcForce(InteractionData idat) {
435 >  void Electrostatic::calcForce(InteractionData &idat) {
436  
437      // utility variables.  Should clean these up and use the Vector3d and
438      // Mat3x3d to replace as many as we can in future versions:
# Line 441 | Line 446 | namespace OpenMD {
446      RealType ct_i, ct_j, ct_ij, a1;
447      RealType riji, ri, ri2, ri3, ri4;
448      RealType pref, vterm, epot, dudr;
449 +    RealType vpair(0.0);
450      RealType scale, sc2;
451      RealType pot_term, preVal, rfVal;
452      RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj;
453      RealType preSw, preSwSc;
454      RealType c1, c2, c3, c4;
455 <    RealType erfcVal, derfcVal;
455 >    RealType erfcVal(1.0), derfcVal(0.0);
456      RealType BigR;
457  
458      Vector3d Q_i, Q_j;
# Line 457 | Line 463 | namespace OpenMD {
463      Vector3d rhatdot2, rhatc4;
464      Vector3d dVdr;
465  
466 +    // variables for indirect (reaction field) interactions for excluded pairs:
467 +    RealType indirect_Pot(0.0);
468 +    RealType indirect_vpair(0.0);
469 +    Vector3d indirect_dVdr(V3Zero);
470 +    Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
471 +
472      pair<RealType, RealType> res;
473      
474      if (!initialized_) initialize();
475      
476 <    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1];
477 <    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2];
476 >    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
477 >    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
478      
479      // some variables we'll need independent of electrostatic type:
480  
481 <    riji = 1.0 / idat.rij;
482 <    Vector3d rhat = idat.d  * riji;
481 >    riji = 1.0 /  *(idat.rij) ;
482 >    Vector3d rhat =  *(idat.d)   * riji;
483  
484      // logicals
485  
# Line 481 | Line 493 | namespace OpenMD {
493      bool j_is_SplitDipole = data2.is_SplitDipole;
494      bool j_is_Quadrupole = data2.is_Quadrupole;
495      
496 <    if (i_is_Charge)
496 >    if (i_is_Charge) {
497        q_i = data1.charge;
498 +      if (idat.excluded) {
499 +        *(idat.skippedCharge2) += q_i;
500 +      }
501 +    }
502  
503      if (i_is_Dipole) {
504        mu_i = data1.dipole_moment;
505 <      uz_i = idat.eFrame1.getColumn(2);
505 >      uz_i = idat.eFrame1->getColumn(2);
506        
507        ct_i = dot(uz_i, rhat);
508  
# Line 502 | Line 518 | namespace OpenMD {
518        qyy_i = Q_i.y();
519        qzz_i = Q_i.z();
520        
521 <      ux_i = idat.eFrame1.getColumn(0);
522 <      uy_i = idat.eFrame1.getColumn(1);
523 <      uz_i = idat.eFrame1.getColumn(2);
521 >      ux_i = idat.eFrame1->getColumn(0);
522 >      uy_i = idat.eFrame1->getColumn(1);
523 >      uz_i = idat.eFrame1->getColumn(2);
524  
525        cx_i = dot(ux_i, rhat);
526        cy_i = dot(uy_i, rhat);
# Line 515 | Line 531 | namespace OpenMD {
531        duduz_i = V3Zero;
532      }
533  
534 <    if (j_is_Charge)
534 >    if (j_is_Charge) {
535        q_j = data2.charge;
536 +      if (idat.excluded) {
537 +        *(idat.skippedCharge1) += q_j;
538 +      }
539 +    }
540  
541 +
542      if (j_is_Dipole) {
543        mu_j = data2.dipole_moment;
544 <      uz_j = idat.eFrame2.getColumn(2);
544 >      uz_j = idat.eFrame2->getColumn(2);
545        
546        ct_j = dot(uz_j, rhat);
547  
# Line 536 | Line 557 | namespace OpenMD {
557        qyy_j = Q_j.y();
558        qzz_j = Q_j.z();
559        
560 <      ux_j = idat.eFrame2.getColumn(0);
561 <      uy_j = idat.eFrame2.getColumn(1);
562 <      uz_j = idat.eFrame2.getColumn(2);
560 >      ux_j = idat.eFrame2->getColumn(0);
561 >      uy_j = idat.eFrame2->getColumn(1);
562 >      uz_j = idat.eFrame2->getColumn(2);
563  
564        cx_j = dot(ux_j, rhat);
565        cy_j = dot(uy_j, rhat);
# Line 557 | Line 578 | namespace OpenMD {
578        if (j_is_Charge) {
579          if (screeningMethod_ == DAMPED) {
580            // assemble the damping variables
581 <          res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
581 >          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
582            erfcVal = res.first;
583            derfcVal = res.second;
584            c1 = erfcVal * riji;
# Line 567 | Line 588 | namespace OpenMD {
588            c2 = c1 * riji;
589          }
590  
591 <        preVal = idat.electroMult * pre11_ * q_i * q_j;
591 >        preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
592          
593          if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
594            vterm = preVal * (c1 - c1c_);
595 <          dudr  = -idat.sw * preVal * c2;
595 >          dudr  = - *(idat.sw)  * preVal * c2;
596  
597          } else if (summationMethod_ == esm_SHIFTED_FORCE)  {
598 <          vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - cutoffRadius_) );
599 <          dudr  = idat.sw * preVal * (c2c_ - c2);
598 >          vterm = preVal * ( c1 - c1c_ + c2c_*( *(idat.rij)  - cutoffRadius_) );
599 >          dudr  =  *(idat.sw)  * preVal * (c2c_ - c2);
600  
601          } else if (summationMethod_ == esm_REACTION_FIELD) {
602 <          rfVal = idat.electroMult * preRF_ * idat.rij * idat.rij;
602 >          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij);
603 >
604            vterm = preVal * ( riji + rfVal );            
605 <          dudr  = idat.sw * preVal * ( 2.0 * rfVal - riji ) * riji;
605 >          dudr  =  *(idat.sw)  * preVal * ( 2.0 * rfVal - riji ) * riji;
606 >          
607 >          // if this is an excluded pair, there are still indirect
608 >          // interactions via the reaction field we must worry about:
609  
610 +          if (idat.excluded) {
611 +            indirect_vpair += preVal * rfVal;
612 +            indirect_Pot += *(idat.sw) * preVal * rfVal;
613 +            indirect_dVdr += *(idat.sw)  * preVal * 2.0 * rfVal  * riji * rhat;
614 +          }
615 +          
616          } else {
586 –          vterm = preVal * riji * erfcVal;            
617  
618 <          dudr  = - idat.sw * preVal * c2;
618 >          vterm = preVal * riji * erfcVal;          
619 >          dudr  = -  *(idat.sw)  * preVal * c2;
620  
621          }
591 –
592 –        idat.vpair += vterm;
593 –        epot += idat.sw * vterm;
622  
623 <        dVdr += dudr * rhat;      
623 >        vpair += vterm;
624 >        epot +=  *(idat.sw)  * vterm;
625 >        dVdr += dudr * rhat;                
626        }
627  
628        if (j_is_Dipole) {
629          // pref is used by all the possible methods
630 <        pref = idat.electroMult * pre12_ * q_i * mu_j;
631 <        preSw = idat.sw * pref;
630 >        pref =  *(idat.electroMult) * pre12_ * q_i * mu_j;
631 >        preSw =  *(idat.sw)  * pref;
632  
633          if (summationMethod_ == esm_REACTION_FIELD) {
634            ri2 = riji * riji;
635            ri3 = ri2 * riji;
636      
637 <          vterm = - pref * ct_j * ( ri2 - preRF2_ * idat.rij );
638 <          idat.vpair += vterm;
639 <          epot += idat.sw * vterm;
637 >          vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
638 >          vpair += vterm;
639 >          epot +=  *(idat.sw)  * vterm;
640  
641            dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
642 <          duduz_j += -preSw * rhat * (ri2 - preRF2_ * idat.rij);  
642 >          duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );  
643  
644 +          // Even if we excluded this pair from direct interactions,
645 +          // we still have the reaction-field-mediated charge-dipole
646 +          // interaction:
647 +
648 +          if (idat.excluded) {
649 +            indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij);
650 +            indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij);
651 +            indirect_dVdr += preSw * preRF2_ * uz_j;
652 +            indirect_duduz_j += preSw * rhat * preRF2_ *  *(idat.rij);
653 +          }
654 +                      
655          } else {
656            // determine the inverse r used if we have split dipoles
657            if (j_is_SplitDipole) {
658 <            BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
658 >            BigR = sqrt( *(idat.r2) + 0.25 * d_j * d_j);
659              ri = 1.0 / BigR;
660 <            scale = idat.rij * ri;
660 >            scale =  *(idat.rij)  * ri;
661            } else {
662              ri = riji;
663              scale = 1.0;
# Line 626 | Line 667 | namespace OpenMD {
667  
668            if (screeningMethod_ == DAMPED) {
669              // assemble the damping variables
670 <            res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
670 >            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
671              erfcVal = res.first;
672              derfcVal = res.second;
673              c1 = erfcVal * ri;
# Line 643 | Line 684 | namespace OpenMD {
684            // calculate the potential
685            pot_term =  scale * c2;
686            vterm = -pref * ct_j * pot_term;
687 <          idat.vpair += vterm;
688 <          epot += idat.sw * vterm;
687 >          vpair += vterm;
688 >          epot +=  *(idat.sw)  * vterm;
689              
690            // calculate derivatives for forces and torques
691  
# Line 659 | Line 700 | namespace OpenMD {
700          cx2 = cx_j * cx_j;
701          cy2 = cy_j * cy_j;
702          cz2 = cz_j * cz_j;
703 <        pref =  idat.electroMult * pre14_ * q_i * one_third_;
703 >        pref =   *(idat.electroMult) * pre14_ * q_i * one_third_;
704            
705          if (screeningMethod_ == DAMPED) {
706            // assemble the damping variables
707 <          res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
707 >          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
708            erfcVal = res.first;
709            derfcVal = res.second;
710            c1 = erfcVal * riji;
# Line 678 | Line 719 | namespace OpenMD {
719          }
720  
721          // precompute variables for convenience
722 <        preSw = idat.sw * pref;
722 >        preSw =  *(idat.sw)  * pref;
723          c2ri = c2 * riji;
724          c3ri = c3 * riji;
725 <        c4rij = c4 * idat.rij;
725 >        c4rij = c4 *  *(idat.rij) ;
726          rhatdot2 = 2.0 * rhat * c3;
727          rhatc4 = rhat * c4rij;
728  
# Line 690 | Line 731 | namespace OpenMD {
731                       qyy_j * (cy2*c3 - c2ri) +
732                       qzz_j * (cz2*c3 - c2ri) );
733          vterm = pref * pot_term;
734 <        idat.vpair += vterm;
735 <        epot += idat.sw * vterm;
734 >        vpair += vterm;
735 >        epot +=  *(idat.sw)  * vterm;
736                  
737          // calculate derivatives for the forces and torques
738  
# Line 709 | Line 750 | namespace OpenMD {
750  
751        if (j_is_Charge) {
752          // variables used by all the methods
753 <        pref = idat.electroMult * pre12_ * q_j * mu_i;
754 <        preSw = idat.sw * pref;
753 >        pref =  *(idat.electroMult) * pre12_ * q_j * mu_i;
754 >        preSw =  *(idat.sw)  * pref;
755  
756          if (summationMethod_ == esm_REACTION_FIELD) {
757  
758            ri2 = riji * riji;
759            ri3 = ri2 * riji;
760  
761 <          vterm = pref * ct_i * ( ri2 - preRF2_ * idat.rij );
762 <          idat.vpair += vterm;
763 <          epot += idat.sw * vterm;
761 >          vterm = pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
762 >          vpair += vterm;
763 >          epot +=  *(idat.sw)  * vterm;
764            
765            dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
766            
767 <          duduz_i += preSw * rhat * (ri2 - preRF2_ * idat.rij);
767 >          duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
768 >
769 >          // Even if we excluded this pair from direct interactions,
770 >          // we still have the reaction-field-mediated charge-dipole
771 >          // interaction:
772 >
773 >          if (idat.excluded) {
774 >            indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij);
775 >            indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij);
776 >            indirect_dVdr += -preSw * preRF2_ * uz_i;
777 >            indirect_duduz_i += -preSw * rhat * preRF2_ *  *(idat.rij);
778 >          }
779              
780          } else {
781            
782            // determine inverse r if we are using split dipoles
783            if (i_is_SplitDipole) {
784 <            BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
784 >            BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i);
785              ri = 1.0 / BigR;
786 <            scale = idat.rij * ri;
786 >            scale =  *(idat.rij)  * ri;
787            } else {
788              ri = riji;
789              scale = 1.0;
# Line 741 | Line 793 | namespace OpenMD {
793              
794            if (screeningMethod_ == DAMPED) {
795              // assemble the damping variables
796 <            res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
796 >            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
797              erfcVal = res.first;
798              derfcVal = res.second;
799              c1 = erfcVal * ri;
# Line 758 | Line 810 | namespace OpenMD {
810            // calculate the potential
811            pot_term = c2 * scale;
812            vterm = pref * ct_i * pot_term;
813 <          idat.vpair += vterm;
814 <          epot += idat.sw * vterm;
813 >          vpair += vterm;
814 >          epot +=  *(idat.sw)  * vterm;
815  
816            // calculate derivatives for the forces and torques
817            dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
# Line 771 | Line 823 | namespace OpenMD {
823          // variables used by all methods
824          ct_ij = dot(uz_i, uz_j);
825  
826 <        pref = idat.electroMult * pre22_ * mu_i * mu_j;
827 <        preSw = idat.sw * pref;
826 >        pref =  *(idat.electroMult) * pre22_ * mu_i * mu_j;
827 >        preSw =  *(idat.sw)  * pref;
828  
829          if (summationMethod_ == esm_REACTION_FIELD) {
830            ri2 = riji * riji;
# Line 781 | Line 833 | namespace OpenMD {
833  
834            vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
835                             preRF2_ * ct_ij );
836 <          idat.vpair += vterm;
837 <          epot += idat.sw * vterm;
836 >          vpair += vterm;
837 >          epot +=  *(idat.sw)  * vterm;
838              
839            a1 = 5.0 * ct_i * ct_j - ct_ij;
840              
# Line 791 | Line 843 | namespace OpenMD {
843            duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
844            duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
845  
846 +          if (idat.excluded) {
847 +            indirect_vpair +=  - pref * preRF2_ * ct_ij;
848 +            indirect_Pot +=    - preSw * preRF2_ * ct_ij;
849 +            indirect_duduz_i += -preSw * preRF2_ * uz_j;
850 +            indirect_duduz_j += -preSw * preRF2_ * uz_i;
851 +          }
852 +
853          } else {
854            
855            if (i_is_SplitDipole) {
856              if (j_is_SplitDipole) {
857 <              BigR = sqrt(idat.r2 + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
857 >              BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
858              } else {
859 <              BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
859 >              BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i);
860              }
861              ri = 1.0 / BigR;
862 <            scale = idat.rij * ri;
862 >            scale =  *(idat.rij)  * ri;
863            } else {
864              if (j_is_SplitDipole) {
865 <              BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
865 >              BigR = sqrt( *(idat.r2) + 0.25 * d_j * d_j);
866                ri = 1.0 / BigR;
867 <              scale = idat.rij * ri;
867 >              scale =  *(idat.rij)  * ri;
868              } else {
869                ri = riji;
870                scale = 1.0;
# Line 813 | Line 872 | namespace OpenMD {
872            }
873            if (screeningMethod_ == DAMPED) {
874              // assemble damping variables
875 <            res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
875 >            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
876              erfcVal = res.first;
877              derfcVal = res.second;
878              c1 = erfcVal * ri;
# Line 835 | Line 894 | namespace OpenMD {
894            preSwSc = preSw * scale;
895            c2ri = c2 * ri;
896            c3ri = c3 * ri;
897 <          c4rij = c4 * idat.rij;
897 >          c4rij = c4 *  *(idat.rij) ;
898  
899            // calculate the potential
900            pot_term = (ct_ij * c2ri - ctidotj * c3);
901            vterm = pref * pot_term;
902 <          idat.vpair += vterm;
903 <          epot += idat.sw * vterm;
902 >          vpair += vterm;
903 >          epot +=  *(idat.sw)  * vterm;
904  
905            // calculate derivatives for the forces and torques
906            dVdr += preSwSc * ( ctidotj * rhat * c4rij  -
# Line 860 | Line 919 | namespace OpenMD {
919          cy2 = cy_i * cy_i;
920          cz2 = cz_i * cz_i;
921  
922 <        pref = idat.electroMult * pre14_ * q_j * one_third_;
922 >        pref =  *(idat.electroMult) * pre14_ * q_j * one_third_;
923  
924          if (screeningMethod_ == DAMPED) {
925            // assemble the damping variables
926 <          res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
926 >          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
927            erfcVal = res.first;
928            derfcVal = res.second;
929            c1 = erfcVal * riji;
# Line 879 | Line 938 | namespace OpenMD {
938          }
939            
940          // precompute some variables for convenience
941 <        preSw = idat.sw * pref;
941 >        preSw =  *(idat.sw)  * pref;
942          c2ri = c2 * riji;
943          c3ri = c3 * riji;
944 <        c4rij = c4 * idat.rij;
944 >        c4rij = c4 *  *(idat.rij) ;
945          rhatdot2 = 2.0 * rhat * c3;
946          rhatc4 = rhat * c4rij;
947  
# Line 892 | Line 951 | namespace OpenMD {
951                       qzz_i * (cz2 * c3 - c2ri) );
952          
953          vterm = pref * pot_term;
954 <        idat.vpair += vterm;
955 <        epot += idat.sw * vterm;
954 >        vpair += vterm;
955 >        epot +=  *(idat.sw)  * vterm;
956  
957          // calculate the derivatives for the forces and torques
958  
# Line 907 | Line 966 | namespace OpenMD {
966        }
967      }
968  
910 –    idat.pot += epot;
911 –    idat.f1 += dVdr;
969  
970 <    if (i_is_Dipole || i_is_Quadrupole)
971 <      idat.t1 -= cross(uz_i, duduz_i);
972 <    if (i_is_Quadrupole) {
973 <      idat.t1 -= cross(ux_i, dudux_i);
917 <      idat.t1 -= cross(uy_i, duduy_i);
918 <    }
919 <
920 <    if (j_is_Dipole || j_is_Quadrupole)
921 <      idat.t2 -= cross(uz_j, duduz_j);
922 <    if (j_is_Quadrupole) {
923 <      idat.t2 -= cross(uz_j, dudux_j);
924 <      idat.t2 -= cross(uz_j, duduy_j);
925 <    }
926 <
927 <    return;
928 <  }  
929 <
930 <  void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) {
931 <
932 <    if (!initialized_) initialize();
933 <    
934 <    ElectrostaticAtomData data1 = ElectrostaticMap[skdat.atype1];
935 <    ElectrostaticAtomData data2 = ElectrostaticMap[skdat.atype2];
936 <    
937 <    // logicals
938 <
939 <    bool i_is_Charge = data1.is_Charge;
940 <    bool i_is_Dipole = data1.is_Dipole;
941 <
942 <    bool j_is_Charge = data2.is_Charge;
943 <    bool j_is_Dipole = data2.is_Dipole;
944 <
945 <    RealType q_i, q_j;
946 <    
947 <    // The skippedCharge computation is needed by the real-space cutoff methods
948 <    // (i.e. shifted force and shifted potential)
949 <
950 <    if (i_is_Charge) {
951 <      q_i = data1.charge;
952 <      skdat.skippedCharge2 += q_i;
953 <    }
954 <
955 <    if (j_is_Charge) {
956 <      q_j = data2.charge;
957 <      skdat.skippedCharge1 += q_j;
958 <    }
959 <
960 <    // the rest of this function should only be necessary for reaction field.
961 <
962 <    if (summationMethod_ == esm_REACTION_FIELD) {
963 <      RealType riji, ri2, ri3;
964 <      RealType q_i, mu_i, ct_i;
965 <      RealType q_j, mu_j, ct_j;
966 <      RealType preVal, rfVal, vterm, dudr, pref, myPot;
967 <      Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
968 <
969 <      // some variables we'll need independent of electrostatic type:
970 >    if (!idat.excluded) {
971 >      *(idat.vpair) += vpair;
972 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
973 >      *(idat.f1) += dVdr;
974        
975 <      riji = 1.0 / skdat.rij;
976 <      rhat = skdat.d  * riji;
977 <
978 <      if (i_is_Dipole) {
979 <        mu_i = data1.dipole_moment;
976 <        uz_i = skdat.eFrame1.getColumn(2);      
977 <        ct_i = dot(uz_i, rhat);
978 <        duduz_i = V3Zero;
975 >      if (i_is_Dipole || i_is_Quadrupole)
976 >        *(idat.t1) -= cross(uz_i, duduz_i);
977 >      if (i_is_Quadrupole) {
978 >        *(idat.t1) -= cross(ux_i, dudux_i);
979 >        *(idat.t1) -= cross(uy_i, duduy_i);
980        }
981 <            
982 <      if (j_is_Dipole) {
983 <        mu_j = data2.dipole_moment;
984 <        uz_j = skdat.eFrame2.getColumn(2);      
985 <        ct_j = dot(uz_j, rhat);
986 <        duduz_j = V3Zero;
981 >      
982 >      if (j_is_Dipole || j_is_Quadrupole)
983 >        *(idat.t2) -= cross(uz_j, duduz_j);
984 >      if (j_is_Quadrupole) {
985 >        *(idat.t2) -= cross(uz_j, dudux_j);
986 >        *(idat.t2) -= cross(uz_j, duduy_j);
987        }
988 <    
989 <      if (i_is_Charge) {
990 <        if (j_is_Charge) {
991 <          preVal = skdat.electroMult * pre11_ * q_i * q_j;
992 <          rfVal = preRF_ * skdat.rij * skdat.rij;
993 <          vterm = preVal * rfVal;
994 <          myPot += skdat.sw * vterm;        
995 <          dudr  = skdat.sw * preVal * 2.0 * rfVal * riji;        
995 <          dVdr += dudr * rhat;
996 <        }
997 <        
998 <        if (j_is_Dipole) {
999 <          ri2 = riji * riji;
1000 <          ri3 = ri2 * riji;        
1001 <          pref = skdat.electroMult * pre12_ * q_i * mu_j;
1002 <          vterm = - pref * ct_j * ( ri2 - preRF2_ * skdat.rij );
1003 <          myPot += skdat.sw * vterm;        
1004 <          dVdr += -skdat.sw * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1005 <          duduz_j += -skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1006 <        }
1007 <      }
1008 <      if (i_is_Dipole) {
1009 <        if (j_is_Charge) {
1010 <          ri2 = riji * riji;
1011 <          ri3 = ri2 * riji;        
1012 <          pref = skdat.electroMult * pre12_ * q_j * mu_i;
1013 <          vterm = - pref * ct_i * ( ri2 - preRF2_ * skdat.rij );
1014 <          myPot += skdat.sw * vterm;        
1015 <          dVdr += skdat.sw * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);      
1016 <          duduz_i += skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1017 <        }
1018 <      }
988 >
989 >    } else {
990 >
991 >      // only accumulate the forces and torques resulting from the
992 >      // indirect reaction field terms.
993 >      *(idat.vpair) += indirect_vpair;
994 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
995 >      *(idat.f1) += indirect_dVdr;
996        
1020 –      // accumulate the forces and torques resulting from the self term
1021 –      skdat.pot += myPot;
1022 –      skdat.f1 += dVdr;
1023 –      
997        if (i_is_Dipole)
998 <        skdat.t1 -= cross(uz_i, duduz_i);
998 >        *(idat.t1) -= cross(uz_i, indirect_duduz_i);
999        if (j_is_Dipole)
1000 <        skdat.t2 -= cross(uz_j, duduz_j);
1000 >        *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1001      }
1002 <  }
1002 >
1003 >
1004 >    return;
1005 >  }  
1006      
1007 <  void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) {
1007 >  void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1008      RealType mu1, preVal, chg1, self;
1009      
1010      if (!initialized_) initialize();
1011 <    
1012 <    ElectrostaticAtomData data = ElectrostaticMap[scdat.atype];
1011 >
1012 >    ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1013    
1014      // logicals
1039 –
1015      bool i_is_Charge = data.is_Charge;
1016      bool i_is_Dipole = data.is_Dipole;
1017  
# Line 1044 | Line 1019 | namespace OpenMD {
1019        if (i_is_Dipole) {
1020          mu1 = data.dipole_moment;          
1021          preVal = pre22_ * preRF2_ * mu1 * mu1;
1022 <        scdat.pot -= 0.5 * preVal;
1022 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal;
1023          
1024          // The self-correction term adds into the reaction field vector
1025 <        Vector3d uz_i = scdat.eFrame.getColumn(2);
1025 >        Vector3d uz_i = sdat.eFrame->getColumn(2);
1026          Vector3d ei = preVal * uz_i;
1027  
1028          // This looks very wrong.  A vector crossed with itself is zero.
1029 <        scdat.t -= cross(uz_i, ei);
1029 >        *(sdat.t) -= cross(uz_i, ei);
1030        }
1031      } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1032        if (i_is_Charge) {        
1033          chg1 = data.charge;
1034          if (screeningMethod_ == DAMPED) {
1035 <          self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1035 >          self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_;
1036          } else {        
1037 <          self = - 0.5 * rcuti_ * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1037 >          self = - 0.5 * rcuti_ * chg1 * (chg1 +  *(sdat.skippedCharge)) * pre11_;
1038          }
1039 <        scdat.pot += self;
1039 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1040        }
1041      }
1042    }
1043  
1044 <  RealType Electrostatic::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
1044 >  RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
1045      // This seems to work moderately well as a default.  There's no
1046      // inherent scale for 1/r interactions that we can standardize.
1047      // 12 angstroms seems to be a reasonably good guess for most

Diff Legend

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