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 1981 by gezelter, Mon Apr 14 18:32:51 2014 UTC vs.
Revision 1993 by gezelter, Tue Apr 29 17:32:31 2014 UTC

# Line 768 | Line 768 | namespace OpenMD {
768      Tb.zero(); // Torque on site b
769      Ea.zero(); // Electric field at site a
770      Eb.zero(); // Electric field at site b
771 +    Pa = 0.0;  // Site potential at site a
772 +    Pb = 0.0;  // Site potential at site b
773      dUdCa = 0.0; // fluctuating charge force at site a
774      dUdCb = 0.0; // fluctuating charge force at site a
775      
# Line 842 | Line 844 | namespace OpenMD {
844        if (idat.excluded) {
845          *(idat.skippedCharge2) += C_a;
846        } else {
847 <        // only do the field if we're not excluded:
847 >        // only do the field and site potentials if we're not excluded:
848          Eb -= C_a *  pre11_ * dv01 * rhat;
849 +        Pb += C_a *  pre11_ * v01;
850        }
851      }
852      
# Line 851 | Line 854 | namespace OpenMD {
854        D_a = *(idat.dipole1);
855        rdDa = dot(rhat, D_a);
856        rxDa = cross(rhat, D_a);
857 <      if (!idat.excluded)
857 >      if (!idat.excluded) {
858          Eb -=  pre12_ * ((dv11-v11or) * rdDa * rhat + v11or * D_a);
859 +        Pb +=  pre12_ * v11 * rdDa;
860 +      }
861 +
862      }
863      
864      if (a_is_Quadrupole) {
# Line 862 | Line 868 | namespace OpenMD {
868        rQa = rhat * Q_a;
869        rdQar = dot(rhat, Qar);
870        rxQar = cross(rhat, Qar);
871 <      if (!idat.excluded)
871 >      if (!idat.excluded) {
872          Eb -= pre14_ * (trQa * rhat * dv21 + 2.0 * Qar * v22or
873                          + rdQar * rhat * (dv22 - 2.0*v22or));
874 +        Pb += pre14_ * (v21 * trQa + v22 * rdQar);
875 +      }
876      }
877      
878      if (b_is_Charge) {
# Line 878 | Line 886 | namespace OpenMD {
886        } else {
887          // only do the field if we're not excluded:
888          Ea += C_b *  pre11_ * dv01 * rhat;
889 +        Pa += C_b *  pre11_ * v01;
890 +
891        }
892      }
893      
# Line 885 | Line 895 | namespace OpenMD {
895        D_b = *(idat.dipole2);
896        rdDb = dot(rhat, D_b);
897        rxDb = cross(rhat, D_b);
898 <      if (!idat.excluded)
898 >      if (!idat.excluded) {
899          Ea += pre12_ * ((dv11-v11or) * rdDb * rhat + v11or * D_b);
900 +        Pa += pre12_ * v11 * rdDb;
901 +      }
902      }
903      
904      if (b_is_Quadrupole) {
# Line 896 | Line 908 | namespace OpenMD {
908        rQb = rhat * Q_b;
909        rdQbr = dot(rhat, Qbr);
910        rxQbr = cross(rhat, Qbr);
911 <      if (!idat.excluded)
911 >      if (!idat.excluded) {
912          Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or
913                          + rdQbr * rhat * (dv22 - 2.0*v22or));
914 +        Pa += pre14_ * (v21 * trQb + v22 * rdQbr);
915 +      }
916      }
917          
918  
# Line 1097 | Line 1111 | namespace OpenMD {
1111      if (idat.doElectricField) {
1112        *(idat.eField1) += Ea * *(idat.electroMult);
1113        *(idat.eField2) += Eb * *(idat.electroMult);
1114 +    }
1115 +
1116 +    if (idat.doSitePotential) {
1117 +      *(idat.sPot1) += Pa * *(idat.electroMult);
1118 +      *(idat.sPot2) += Pb * *(idat.electroMult);
1119      }
1120  
1121      if (a_is_Fluctuating) *(idat.dVdFQ1) += dUdCa * *(idat.sw);

Diff Legend

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