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 1584 by gezelter, Fri Jun 17 20:16:35 2011 UTC

# 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),
# Line 56 | Line 56 | namespace OpenMD {
56    
57    void Electrostatic::initialize() {
58  
59 +    Globals* simParams_ = info_->getSimParams();
60 +
61      summationMap_["HARD"]               = esm_HARD;
62      summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
63      summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
# Line 111 | Line 113 | namespace OpenMD {
113        } else {
114          // throw error
115          sprintf( painCave.errMsg,
116 <                 "SimInfo error: Unknown electrostaticSummationMethod.\n"
116 >                 "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
117                   "\t(Input file specified %s .)\n"
118                   "\telectrostaticSummationMethod must be one of: \"none\",\n"
119                   "\t\"shifted_potential\", \"shifted_force\", or \n"
# Line 405 | Line 407 | namespace OpenMD {
407      return;
408    }
409    
410 <  void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
411 <                                                    RealType theRSW ) {
410 <    cutoffRadius_ = theECR;
410 >  void Electrostatic::setCutoffRadius( RealType rCut ) {
411 >    cutoffRadius_ = rCut;
412      rrf_ = cutoffRadius_;
412 –    rt_ = theRSW;
413      haveCutoffRadius_ = true;
414    }
415 +
416 +  void Electrostatic::setSwitchingRadius( RealType rSwitch ) {
417 +    rt_ = rSwitch;
418 +  }
419    void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
420      summationMethod_ = esm;
421    }
# Line 427 | Line 431 | namespace OpenMD {
431      haveDielectric_ = true;
432    }
433  
434 <  void Electrostatic::calcForce(InteractionData idat) {
434 >  void Electrostatic::calcForce(InteractionData &idat) {
435  
436      // utility variables.  Should clean these up and use the Vector3d and
437      // Mat3x3d to replace as many as we can in future versions:
# Line 461 | Line 465 | namespace OpenMD {
465      
466      if (!initialized_) initialize();
467      
468 <    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1];
469 <    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2];
468 >    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
469 >    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
470      
471      // some variables we'll need independent of electrostatic type:
472  
473 <    riji = 1.0 / idat.rij;
474 <    Vector3d rhat = idat.d  * riji;
473 >    riji = 1.0 /  *(idat.rij) ;
474 >    Vector3d rhat =  *(idat.d)   * riji;
475  
476      // logicals
477  
# Line 486 | Line 490 | namespace OpenMD {
490  
491      if (i_is_Dipole) {
492        mu_i = data1.dipole_moment;
493 <      uz_i = idat.eFrame1.getColumn(2);
493 >      uz_i = idat.eFrame1->getColumn(2);
494        
495        ct_i = dot(uz_i, rhat);
496  
# Line 502 | Line 506 | namespace OpenMD {
506        qyy_i = Q_i.y();
507        qzz_i = Q_i.z();
508        
509 <      ux_i = idat.eFrame1.getColumn(0);
510 <      uy_i = idat.eFrame1.getColumn(1);
511 <      uz_i = idat.eFrame1.getColumn(2);
509 >      ux_i = idat.eFrame1->getColumn(0);
510 >      uy_i = idat.eFrame1->getColumn(1);
511 >      uz_i = idat.eFrame1->getColumn(2);
512  
513        cx_i = dot(ux_i, rhat);
514        cy_i = dot(uy_i, rhat);
# Line 520 | Line 524 | namespace OpenMD {
524  
525      if (j_is_Dipole) {
526        mu_j = data2.dipole_moment;
527 <      uz_j = idat.eFrame2.getColumn(2);
527 >      uz_j = idat.eFrame2->getColumn(2);
528        
529        ct_j = dot(uz_j, rhat);
530  
# Line 536 | Line 540 | namespace OpenMD {
540        qyy_j = Q_j.y();
541        qzz_j = Q_j.z();
542        
543 <      ux_j = idat.eFrame2.getColumn(0);
544 <      uy_j = idat.eFrame2.getColumn(1);
545 <      uz_j = idat.eFrame2.getColumn(2);
543 >      ux_j = idat.eFrame2->getColumn(0);
544 >      uy_j = idat.eFrame2->getColumn(1);
545 >      uz_j = idat.eFrame2->getColumn(2);
546  
547        cx_j = dot(ux_j, rhat);
548        cy_j = dot(uy_j, rhat);
# Line 557 | Line 561 | namespace OpenMD {
561        if (j_is_Charge) {
562          if (screeningMethod_ == DAMPED) {
563            // assemble the damping variables
564 <          res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
564 >          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
565            erfcVal = res.first;
566            derfcVal = res.second;
567            c1 = erfcVal * riji;
# Line 567 | Line 571 | namespace OpenMD {
571            c2 = c1 * riji;
572          }
573  
574 <        preVal = idat.electroMult * pre11_ * q_i * q_j;
574 >        preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
575          
576          if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
577            vterm = preVal * (c1 - c1c_);
578 <          dudr  = -idat.sw * preVal * c2;
578 >          dudr  = - *(idat.sw)  * preVal * c2;
579  
580          } else if (summationMethod_ == esm_SHIFTED_FORCE)  {
581 <          vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - cutoffRadius_) );
582 <          dudr  = idat.sw * preVal * (c2c_ - c2);
581 >          vterm = preVal * ( c1 - c1c_ + c2c_*( *(idat.rij)  - cutoffRadius_) );
582 >          dudr  =  *(idat.sw)  * preVal * (c2c_ - c2);
583  
584          } else if (summationMethod_ == esm_REACTION_FIELD) {
585 <          rfVal = idat.electroMult * preRF_ * idat.rij * idat.rij;
585 >          rfVal =  *(idat.electroMult) * preRF_ *  *(idat.rij)  *  *(idat.rij) ;
586            vterm = preVal * ( riji + rfVal );            
587 <          dudr  = idat.sw * preVal * ( 2.0 * rfVal - riji ) * riji;
587 >          dudr  =  *(idat.sw)  * preVal * ( 2.0 * rfVal - riji ) * riji;
588  
589          } else {
590            vterm = preVal * riji * erfcVal;            
591  
592 <          dudr  = - idat.sw * preVal * c2;
592 >          dudr  = -  *(idat.sw)  * preVal * c2;
593  
594          }
595  
596 <        idat.vpair += vterm;
597 <        epot += idat.sw * vterm;
596 >        *(idat.vpair) += vterm;
597 >        epot +=  *(idat.sw)  * vterm;
598  
599          dVdr += dudr * rhat;      
600        }
601  
602        if (j_is_Dipole) {
603          // pref is used by all the possible methods
604 <        pref = idat.electroMult * pre12_ * q_i * mu_j;
605 <        preSw = idat.sw * pref;
604 >        pref =  *(idat.electroMult) * pre12_ * q_i * mu_j;
605 >        preSw =  *(idat.sw)  * pref;
606  
607          if (summationMethod_ == esm_REACTION_FIELD) {
608            ri2 = riji * riji;
609            ri3 = ri2 * riji;
610      
611 <          vterm = - pref * ct_j * ( ri2 - preRF2_ * idat.rij );
612 <          idat.vpair += vterm;
613 <          epot += idat.sw * vterm;
611 >          vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
612 >          *(idat.vpair) += vterm;
613 >          epot +=  *(idat.sw)  * vterm;
614  
615            dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
616 <          duduz_j += -preSw * rhat * (ri2 - preRF2_ * idat.rij);  
616 >          duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );  
617  
618          } else {
619            // determine the inverse r used if we have split dipoles
620            if (j_is_SplitDipole) {
621 <            BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
621 >            BigR = sqrt( *(idat.r2) + 0.25 * d_j * d_j);
622              ri = 1.0 / BigR;
623 <            scale = idat.rij * ri;
623 >            scale =  *(idat.rij)  * ri;
624            } else {
625              ri = riji;
626              scale = 1.0;
# Line 626 | Line 630 | namespace OpenMD {
630  
631            if (screeningMethod_ == DAMPED) {
632              // assemble the damping variables
633 <            res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
633 >            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
634              erfcVal = res.first;
635              derfcVal = res.second;
636              c1 = erfcVal * ri;
# Line 643 | Line 647 | namespace OpenMD {
647            // calculate the potential
648            pot_term =  scale * c2;
649            vterm = -pref * ct_j * pot_term;
650 <          idat.vpair += vterm;
651 <          epot += idat.sw * vterm;
650 >          *(idat.vpair) += vterm;
651 >          epot +=  *(idat.sw)  * vterm;
652              
653            // calculate derivatives for forces and torques
654  
# Line 659 | Line 663 | namespace OpenMD {
663          cx2 = cx_j * cx_j;
664          cy2 = cy_j * cy_j;
665          cz2 = cz_j * cz_j;
666 <        pref =  idat.electroMult * pre14_ * q_i * one_third_;
666 >        pref =   *(idat.electroMult) * pre14_ * q_i * one_third_;
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 * riji;
# Line 678 | Line 682 | namespace OpenMD {
682          }
683  
684          // precompute variables for convenience
685 <        preSw = idat.sw * pref;
685 >        preSw =  *(idat.sw)  * pref;
686          c2ri = c2 * riji;
687          c3ri = c3 * riji;
688 <        c4rij = c4 * idat.rij;
688 >        c4rij = c4 *  *(idat.rij) ;
689          rhatdot2 = 2.0 * rhat * c3;
690          rhatc4 = rhat * c4rij;
691  
# Line 690 | Line 694 | namespace OpenMD {
694                       qyy_j * (cy2*c3 - c2ri) +
695                       qzz_j * (cz2*c3 - c2ri) );
696          vterm = pref * pot_term;
697 <        idat.vpair += vterm;
698 <        epot += idat.sw * vterm;
697 >        *(idat.vpair) += vterm;
698 >        epot +=  *(idat.sw)  * vterm;
699                  
700          // calculate derivatives for the forces and torques
701  
# Line 709 | Line 713 | namespace OpenMD {
713  
714        if (j_is_Charge) {
715          // variables used by all the methods
716 <        pref = idat.electroMult * pre12_ * q_j * mu_i;
717 <        preSw = idat.sw * pref;
716 >        pref =  *(idat.electroMult) * pre12_ * q_j * mu_i;
717 >        preSw =  *(idat.sw)  * pref;
718  
719          if (summationMethod_ == esm_REACTION_FIELD) {
720  
721            ri2 = riji * riji;
722            ri3 = ri2 * riji;
723  
724 <          vterm = pref * ct_i * ( ri2 - preRF2_ * idat.rij );
725 <          idat.vpair += vterm;
726 <          epot += idat.sw * vterm;
724 >          vterm = pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
725 >          *(idat.vpair) += vterm;
726 >          epot +=  *(idat.sw)  * vterm;
727            
728            dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
729            
730 <          duduz_i += preSw * rhat * (ri2 - preRF2_ * idat.rij);
730 >          duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
731              
732          } else {
733            
734            // determine inverse r if we are using split dipoles
735            if (i_is_SplitDipole) {
736 <            BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
736 >            BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i);
737              ri = 1.0 / BigR;
738 <            scale = idat.rij * ri;
738 >            scale =  *(idat.rij)  * ri;
739            } else {
740              ri = riji;
741              scale = 1.0;
# Line 741 | Line 745 | namespace OpenMD {
745              
746            if (screeningMethod_ == DAMPED) {
747              // assemble the damping variables
748 <            res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
748 >            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
749              erfcVal = res.first;
750              derfcVal = res.second;
751              c1 = erfcVal * ri;
# Line 758 | Line 762 | namespace OpenMD {
762            // calculate the potential
763            pot_term = c2 * scale;
764            vterm = pref * ct_i * pot_term;
765 <          idat.vpair += vterm;
766 <          epot += idat.sw * vterm;
765 >          *(idat.vpair) += vterm;
766 >          epot +=  *(idat.sw)  * vterm;
767  
768            // calculate derivatives for the forces and torques
769            dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
# Line 771 | Line 775 | namespace OpenMD {
775          // variables used by all methods
776          ct_ij = dot(uz_i, uz_j);
777  
778 <        pref = idat.electroMult * pre22_ * mu_i * mu_j;
779 <        preSw = idat.sw * pref;
778 >        pref =  *(idat.electroMult) * pre22_ * mu_i * mu_j;
779 >        preSw =  *(idat.sw)  * pref;
780  
781          if (summationMethod_ == esm_REACTION_FIELD) {
782            ri2 = riji * riji;
# Line 781 | Line 785 | namespace OpenMD {
785  
786            vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
787                             preRF2_ * ct_ij );
788 <          idat.vpair += vterm;
789 <          epot += idat.sw * vterm;
788 >          *(idat.vpair) += vterm;
789 >          epot +=  *(idat.sw)  * vterm;
790              
791            a1 = 5.0 * ct_i * ct_j - ct_ij;
792              
# Line 795 | Line 799 | namespace OpenMD {
799            
800            if (i_is_SplitDipole) {
801              if (j_is_SplitDipole) {
802 <              BigR = sqrt(idat.r2 + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
802 >              BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
803              } else {
804 <              BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
804 >              BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i);
805              }
806              ri = 1.0 / BigR;
807 <            scale = idat.rij * ri;
807 >            scale =  *(idat.rij)  * ri;
808            } else {
809              if (j_is_SplitDipole) {
810 <              BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
810 >              BigR = sqrt( *(idat.r2) + 0.25 * d_j * d_j);
811                ri = 1.0 / BigR;
812 <              scale = idat.rij * ri;
812 >              scale =  *(idat.rij)  * ri;
813              } else {
814                ri = riji;
815                scale = 1.0;
# Line 813 | Line 817 | namespace OpenMD {
817            }
818            if (screeningMethod_ == DAMPED) {
819              // assemble damping variables
820 <            res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
820 >            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
821              erfcVal = res.first;
822              derfcVal = res.second;
823              c1 = erfcVal * ri;
# Line 835 | Line 839 | namespace OpenMD {
839            preSwSc = preSw * scale;
840            c2ri = c2 * ri;
841            c3ri = c3 * ri;
842 <          c4rij = c4 * idat.rij;
842 >          c4rij = c4 *  *(idat.rij) ;
843  
844            // calculate the potential
845            pot_term = (ct_ij * c2ri - ctidotj * c3);
846            vterm = pref * pot_term;
847 <          idat.vpair += vterm;
848 <          epot += idat.sw * vterm;
847 >          *(idat.vpair) += vterm;
848 >          epot +=  *(idat.sw)  * vterm;
849  
850            // calculate derivatives for the forces and torques
851            dVdr += preSwSc * ( ctidotj * rhat * c4rij  -
# Line 860 | Line 864 | namespace OpenMD {
864          cy2 = cy_i * cy_i;
865          cz2 = cz_i * cz_i;
866  
867 <        pref = idat.electroMult * pre14_ * q_j * one_third_;
867 >        pref =  *(idat.electroMult) * pre14_ * q_j * one_third_;
868  
869          if (screeningMethod_ == DAMPED) {
870            // assemble the damping variables
871 <          res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
871 >          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
872            erfcVal = res.first;
873            derfcVal = res.second;
874            c1 = erfcVal * riji;
# Line 879 | Line 883 | namespace OpenMD {
883          }
884            
885          // precompute some variables for convenience
886 <        preSw = idat.sw * pref;
886 >        preSw =  *(idat.sw)  * pref;
887          c2ri = c2 * riji;
888          c3ri = c3 * riji;
889 <        c4rij = c4 * idat.rij;
889 >        c4rij = c4 *  *(idat.rij) ;
890          rhatdot2 = 2.0 * rhat * c3;
891          rhatc4 = rhat * c4rij;
892  
# Line 892 | Line 896 | namespace OpenMD {
896                       qzz_i * (cz2 * c3 - c2ri) );
897          
898          vterm = pref * pot_term;
899 <        idat.vpair += vterm;
900 <        epot += idat.sw * vterm;
899 >        *(idat.vpair) += vterm;
900 >        epot +=  *(idat.sw)  * vterm;
901  
902          // calculate the derivatives for the forces and torques
903  
# Line 907 | Line 911 | namespace OpenMD {
911        }
912      }
913  
914 <    idat.pot += epot;
915 <    idat.f1 += dVdr;
914 >    (*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
915 >    *(idat.f1) += dVdr;
916  
917      if (i_is_Dipole || i_is_Quadrupole)
918 <      idat.t1 -= cross(uz_i, duduz_i);
918 >      *(idat.t1) -= cross(uz_i, duduz_i);
919      if (i_is_Quadrupole) {
920 <      idat.t1 -= cross(ux_i, dudux_i);
921 <      idat.t1 -= cross(uy_i, duduy_i);
920 >      *(idat.t1) -= cross(ux_i, dudux_i);
921 >      *(idat.t1) -= cross(uy_i, duduy_i);
922      }
923 <
923 >    
924      if (j_is_Dipole || j_is_Quadrupole)
925 <      idat.t2 -= cross(uz_j, duduz_j);
925 >      *(idat.t2) -= cross(uz_j, duduz_j);
926      if (j_is_Quadrupole) {
927 <      idat.t2 -= cross(uz_j, dudux_j);
928 <      idat.t2 -= cross(uz_j, duduy_j);
927 >      *(idat.t2) -= cross(uz_j, dudux_j);
928 >      *(idat.t2) -= cross(uz_j, duduy_j);
929      }
930  
931      return;
932    }  
933  
934 <  void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) {
934 >  void Electrostatic::calcSkipCorrection(InteractionData &idat) {
935  
936      if (!initialized_) initialize();
937      
938 <    ElectrostaticAtomData data1 = ElectrostaticMap[skdat.atype1];
939 <    ElectrostaticAtomData data2 = ElectrostaticMap[skdat.atype2];
938 >    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
939 >    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
940      
941      // logicals
942  
# Line 949 | Line 953 | namespace OpenMD {
953  
954      if (i_is_Charge) {
955        q_i = data1.charge;
956 <      skdat.skippedCharge2 += q_i;
956 >      *(idat.skippedCharge2) += q_i;
957      }
958  
959      if (j_is_Charge) {
960        q_j = data2.charge;
961 <      skdat.skippedCharge1 += q_j;
961 >      *(idat.skippedCharge1) += q_j;
962      }
963  
964      // the rest of this function should only be necessary for reaction field.
965  
966      if (summationMethod_ == esm_REACTION_FIELD) {
967        RealType riji, ri2, ri3;
968 <      RealType q_i, mu_i, ct_i;
969 <      RealType q_j, mu_j, ct_j;
970 <      RealType preVal, rfVal, vterm, dudr, pref, myPot;
968 >      RealType mu_i, ct_i;
969 >      RealType mu_j, ct_j;
970 >      RealType preVal, rfVal, vterm, dudr, pref, myPot(0.0);
971        Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
972  
973        // some variables we'll need independent of electrostatic type:
974        
975 <      riji = 1.0 / skdat.rij;
976 <      rhat = skdat.d  * riji;
975 >      riji = 1.0 /  *(idat.rij) ;
976 >      rhat =  *(idat.d)  * riji;
977  
978        if (i_is_Dipole) {
979          mu_i = data1.dipole_moment;
980 <        uz_i = skdat.eFrame1.getColumn(2);      
980 >        uz_i = idat.eFrame1->getColumn(2);      
981          ct_i = dot(uz_i, rhat);
982          duduz_i = V3Zero;
983        }
984              
985        if (j_is_Dipole) {
986          mu_j = data2.dipole_moment;
987 <        uz_j = skdat.eFrame2.getColumn(2);      
987 >        uz_j = idat.eFrame2->getColumn(2);      
988          ct_j = dot(uz_j, rhat);
989          duduz_j = V3Zero;
990        }
991      
992        if (i_is_Charge) {
993          if (j_is_Charge) {
994 <          preVal = skdat.electroMult * pre11_ * q_i * q_j;
995 <          rfVal = preRF_ * skdat.rij * skdat.rij;
994 >          preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
995 >          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij) ;
996            vterm = preVal * rfVal;
997 <          myPot += skdat.sw * vterm;        
998 <          dudr  = skdat.sw * preVal * 2.0 * rfVal * riji;        
997 >          myPot +=  *(idat.sw)  * vterm;        
998 >          dudr  =  *(idat.sw)  * preVal * 2.0 * rfVal * riji;        
999            dVdr += dudr * rhat;
1000          }
1001          
1002          if (j_is_Dipole) {
1003            ri2 = riji * riji;
1004            ri3 = ri2 * riji;        
1005 <          pref = skdat.electroMult * pre12_ * q_i * mu_j;
1006 <          vterm = - pref * ct_j * ( ri2 - preRF2_ * skdat.rij );
1007 <          myPot += skdat.sw * vterm;        
1008 <          dVdr += -skdat.sw * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1009 <          duduz_j += -skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1005 >          pref =  *(idat.electroMult) * pre12_ * q_i * mu_j;
1006 >          vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
1007 >          myPot +=  *(idat.sw)  * vterm;        
1008 >          dVdr += - *(idat.sw)  * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1009 >          duduz_j += - *(idat.sw)  * pref * rhat * (ri2 - preRF2_ *  *(idat.rij) );
1010          }
1011        }
1012        if (i_is_Dipole) {
1013          if (j_is_Charge) {
1014            ri2 = riji * riji;
1015            ri3 = ri2 * riji;        
1016 <          pref = skdat.electroMult * pre12_ * q_j * mu_i;
1017 <          vterm = - pref * ct_i * ( ri2 - preRF2_ * skdat.rij );
1018 <          myPot += skdat.sw * vterm;        
1019 <          dVdr += skdat.sw * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);      
1020 <          duduz_i += skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1016 >          pref =  *(idat.electroMult) * pre12_ * q_j * mu_i;
1017 >          vterm = - pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
1018 >          myPot +=  *(idat.sw)  * vterm;        
1019 >          dVdr +=  *(idat.sw)  * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);      
1020 >          duduz_i +=  *(idat.sw)  * pref * rhat * (ri2 - preRF2_ *  *(idat.rij));
1021          }
1022        }
1023        
1024        // accumulate the forces and torques resulting from the self term
1025 <      skdat.pot += myPot;
1026 <      skdat.f1 += dVdr;
1025 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += myPot;
1026 >      *(idat.f1) += dVdr;
1027        
1028        if (i_is_Dipole)
1029 <        skdat.t1 -= cross(uz_i, duduz_i);
1029 >        *(idat.t1) -= cross(uz_i, duduz_i);
1030        if (j_is_Dipole)
1031 <        skdat.t2 -= cross(uz_j, duduz_j);
1031 >        *(idat.t2) -= cross(uz_j, duduz_j);
1032      }
1033    }
1034      
1035 <  void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) {
1035 >  void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1036      RealType mu1, preVal, chg1, self;
1037      
1038      if (!initialized_) initialize();
1039      
1040 <    ElectrostaticAtomData data = ElectrostaticMap[scdat.atype];
1040 >    ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1041    
1042      // logicals
1043  
# Line 1044 | Line 1048 | namespace OpenMD {
1048        if (i_is_Dipole) {
1049          mu1 = data.dipole_moment;          
1050          preVal = pre22_ * preRF2_ * mu1 * mu1;
1051 <        scdat.pot -= 0.5 * preVal;
1051 >        sdat.pot[2] -= 0.5 * preVal;
1052          
1053          // The self-correction term adds into the reaction field vector
1054 <        Vector3d uz_i = scdat.eFrame.getColumn(2);
1054 >        Vector3d uz_i = sdat.eFrame->getColumn(2);
1055          Vector3d ei = preVal * uz_i;
1056  
1057          // This looks very wrong.  A vector crossed with itself is zero.
1058 <        scdat.t -= cross(uz_i, ei);
1058 >        *(sdat.t) -= cross(uz_i, ei);
1059        }
1060      } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1061        if (i_is_Charge) {        
1062          chg1 = data.charge;
1063          if (screeningMethod_ == DAMPED) {
1064 <          self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1064 >          self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_;
1065          } else {        
1066 <          self = - 0.5 * rcuti_ * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1066 >          self = - 0.5 * rcuti_ * chg1 * (chg1 +  *(sdat.skippedCharge)) * pre11_;
1067          }
1068 <        scdat.pot += self;
1068 >        sdat.pot[ELECTROSTATIC_FAMILY] += self;
1069        }
1070      }
1071    }
1072  
1073 <  RealType Electrostatic::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
1073 >  RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
1074      // This seems to work moderately well as a default.  There's no
1075      // inherent scale for 1/r interactions that we can standardize.
1076      // 12 angstroms seems to be a reasonably good guess for most

Diff Legend

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