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 1925 by gezelter, Wed Aug 7 15:24:16 2013 UTC vs.
Revision 1933 by gezelter, Fri Aug 23 15:59:23 2013 UTC

# Line 765 | Line 765 | namespace OpenMD {
765  
766      // Excluded potential that is still computed for fluctuating charges
767      excluded_Pot= 0.0;
768
768  
769      // some variables we'll need independent of electrostatic type:
770  
# Line 887 | Line 886 | namespace OpenMD {
886          Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or
887                          + rdQbr * rhat * (dv22 - 2.0*v22or));
888      }
889 <    
889 >        
890 >
891      if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
892        J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
893      }    
894 <    
894 >
895      if (a_is_Charge) {    
896        
897        if (b_is_Charge) {
898          pref =  pre11_ * *(idat.electroMult);      
899          U  += C_a * C_b * pref * v01;
900          F  += C_a * C_b * pref * dv01 * rhat;
901 <        
901 >
902          // If this is an excluded pair, there are still indirect
903          // interactions via the reaction field we must worry about:
904  
# Line 907 | Line 907 | namespace OpenMD {
907            indirect_Pot += rfContrib;
908            indirect_F   += rfContrib * 2.0 * ri * rhat;
909          }
910 <        
910 >
911          // Fluctuating charge forces are handled via Coulomb integrals
912          // for excluded pairs (i.e. those connected via bonds) and
913          // with the standard charge-charge interaction otherwise.
914  
915 <        if (idat.excluded) {          
915 >        if (idat.excluded) {
916            if (a_is_Fluctuating || b_is_Fluctuating) {
917              coulInt = J->getValueAt( *(idat.rij) );
918 <            if (a_is_Fluctuating)  dUdCa += coulInt * C_b;
919 <            if (b_is_Fluctuating)  dUdCb += coulInt * C_a;
920 <            excluded_Pot += C_a * C_b * coulInt;
921 <          }          
918 >            if (a_is_Fluctuating) dUdCa += C_b * coulInt;
919 >            if (b_is_Fluctuating) dUdCb += C_a * coulInt;          
920 >          }
921          } else {
922            if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
923 <          if (a_is_Fluctuating) dUdCb += C_a * pref * v01;
924 <        }
923 >          if (b_is_Fluctuating) dUdCb += C_a * pref * v01;
924 >        }              
925        }
926  
927        if (b_is_Dipole) {
# Line 988 | Line 987 | namespace OpenMD {
987          F  -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat;
988          Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
989          Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
991
990          // Even if we excluded this pair from direct interactions, we
991          // still have the reaction-field-mediated dipole-dipole
992          // interaction:
# Line 1048 | Line 1046 | namespace OpenMD {
1046          trQaQb = QaQb.trace();
1047          rQaQb = rhat * QaQb;
1048          QaQbr = QaQb * rhat;
1049 <        QaxQb = cross(Q_a, Q_b);
1049 >        QaxQb = mCross(Q_a, Q_b);
1050          rQaQbr = dot(rQa, Qbr);
1051          rQaxQbr = cross(rQa, Qbr);
1052          
# Line 1079 | Line 1077 | namespace OpenMD {
1077          //             + 4.0 * cross(rhat, QbQar)
1078  
1079          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1082
1080        }
1081      }
1082  
# Line 1142 | Line 1139 | namespace OpenMD {
1139          
1140      if (i_is_Fluctuating) {
1141        C_a += *(sdat.flucQ);
1142 <      // dVdFQ is really a force, so this is negative the derivative
1146 <      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1142 >      *(sdat.flucQfrc) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1143        (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1144          (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1145      }
# Line 1441 | Line 1437 | namespace OpenMD {
1437                    dks[i] = dk * skr[i];
1438                  }
1439                  if (data.is_Quadrupole) {
1440 <                  Q = atom->getQuadrupole();
1441 <                  Q *= mPoleConverter;
1446 <                  Qk = Q * kVec;
1440 >                  Q = atom->getQuadrupole() * mPoleConverter;
1441 >                  Qk = Q * kVec;                  
1442                    qk = dot(kVec, Qk);
1443                    qxk[i] = cross(kVec, Qk);
1444                    qkc[i] = qk * ckr[i];
# Line 1504 | Line 1499 | namespace OpenMD {
1499                  RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
1500                                           -ckr[i]*(ckss+dkcs-qkss));
1501                  RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
1502 <                                             +skr[i]*(ckss+dkcs-qkss));
1502 >                                            +skr[i]*(ckss+dkcs-qkss));
1503                
1504                  atom->addFrc( 4.0 * rvol * qfrc * kVec );
1505                  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines