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 1924 by gezelter, Mon Aug 5 21:46:11 2013 UTC vs.
Revision 1981 by gezelter, Mon Apr 14 18:32:51 2014 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  
# Line 57 | Line 61
61   #include "math/erfc.hpp"
62   #include "math/SquareMatrix.hpp"
63   #include "primitives/Molecule.hpp"
64 + #include "flucq/FluctuatingChargeForces.hpp"
65  
61 –
66   namespace OpenMD {
67    
68    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
# Line 67 | Line 71 | namespace OpenMD {
71                                    haveDampingAlpha_(false),
72                                    haveDielectric_(false),
73                                    haveElectroSplines_(false)
74 <  {}
74 >  {
75 >    flucQ_ = new FluctuatingChargeForces(info_);
76 >  }
77    
78 +  void Electrostatic::setForceField(ForceField *ff) {
79 +    forceField_ = ff;
80 +    flucQ_->setForceField(forceField_);
81 +  }
82 +
83 +  void Electrostatic::setSimulatedAtomTypes(set<AtomType*> &simtypes) {
84 +    simTypes_ = simtypes;
85 +    flucQ_->setSimulatedAtomTypes(simTypes_);
86 +  }
87 +
88    void Electrostatic::initialize() {
89      
90      Globals* simParams_ = info_->getSimParams();
# Line 763 | Line 779 | namespace OpenMD {
779  
780      // Excluded potential that is still computed for fluctuating charges
781      excluded_Pot= 0.0;
766 –
782  
783      // some variables we'll need independent of electrostatic type:
784  
# Line 885 | Line 900 | namespace OpenMD {
900          Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or
901                          + rdQbr * rhat * (dv22 - 2.0*v22or));
902      }
903 <    
903 >        
904 >
905      if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
906        J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
907      }    
908 <    
908 >
909      if (a_is_Charge) {    
910        
911        if (b_is_Charge) {
912          pref =  pre11_ * *(idat.electroMult);      
913          U  += C_a * C_b * pref * v01;
914          F  += C_a * C_b * pref * dv01 * rhat;
915 <        
915 >
916          // If this is an excluded pair, there are still indirect
917          // interactions via the reaction field we must worry about:
918  
# Line 905 | Line 921 | namespace OpenMD {
921            indirect_Pot += rfContrib;
922            indirect_F   += rfContrib * 2.0 * ri * rhat;
923          }
924 <        
924 >
925          // Fluctuating charge forces are handled via Coulomb integrals
926          // for excluded pairs (i.e. those connected via bonds) and
927          // with the standard charge-charge interaction otherwise.
928  
929 <        if (idat.excluded) {          
929 >        if (idat.excluded) {
930            if (a_is_Fluctuating || b_is_Fluctuating) {
931              coulInt = J->getValueAt( *(idat.rij) );
932 <            if (a_is_Fluctuating)  dUdCa += coulInt * C_b;
933 <            if (b_is_Fluctuating)  dUdCb += coulInt * C_a;
934 <            excluded_Pot += C_a * C_b * coulInt;
919 <          }          
932 >            if (a_is_Fluctuating) dUdCa += C_b * coulInt;
933 >            if (b_is_Fluctuating) dUdCb += C_a * coulInt;          
934 >          }
935          } else {
936            if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
937 <          if (a_is_Fluctuating) dUdCb += C_a * pref * v01;
938 <        }
937 >          if (b_is_Fluctuating) dUdCb += C_a * pref * v01;
938 >        }              
939        }
940  
941        if (b_is_Dipole) {
# Line 986 | Line 1001 | namespace OpenMD {
1001          F  -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat;
1002          Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
1003          Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
989 –
1004          // Even if we excluded this pair from direct interactions, we
1005          // still have the reaction-field-mediated dipole-dipole
1006          // interaction:
# Line 1046 | Line 1060 | namespace OpenMD {
1060          trQaQb = QaQb.trace();
1061          rQaQb = rhat * QaQb;
1062          QaQbr = QaQb * rhat;
1063 <        QaxQb = cross(Q_a, Q_b);
1063 >        QaxQb = mCross(Q_a, Q_b);
1064          rQaQbr = dot(rQa, Qbr);
1065          rQaxQbr = cross(rQa, Qbr);
1066          
# Line 1077 | Line 1091 | namespace OpenMD {
1091          //             + 4.0 * cross(rhat, QbQar)
1092  
1093          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1080 –
1094        }
1095      }
1096  
# Line 1140 | Line 1153 | namespace OpenMD {
1153          
1154      if (i_is_Fluctuating) {
1155        C_a += *(sdat.flucQ);
1156 <      // dVdFQ is really a force, so this is negative the derivative
1157 <      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1158 <      (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1159 <        (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1156 >
1157 >      flucQ_->getSelfInteraction(sdat.atid, *(sdat.flucQ),  
1158 >                                 (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY],
1159 >                                 *(sdat.flucQfrc) );
1160 >
1161      }
1162  
1163      switch (summationMethod_) {
# Line 1194 | Line 1208 | namespace OpenMD {
1208    }
1209  
1210  
1211 <  void Electrostatic::ReciprocalSpaceSum(potVec& pot) {
1211 >  void Electrostatic::ReciprocalSpaceSum(RealType& pot) {
1212      
1213      RealType kPot = 0.0;
1214      RealType kVir = 0.0;
# Line 1240 | Line 1254 | namespace OpenMD {
1254      
1255      // Calculate and store exponential factors  
1256      
1257 <    vector<vector<Vector3d> > eCos;
1258 <    vector<vector<Vector3d> > eSin;
1257 >    vector<vector<RealType> > elc;
1258 >    vector<vector<RealType> > emc;
1259 >    vector<vector<RealType> > enc;
1260 >    vector<vector<RealType> > els;
1261 >    vector<vector<RealType> > ems;
1262 >    vector<vector<RealType> > ens;
1263      
1264      int nMax = info_->getNAtoms();
1265      
1266 <    eCos.resize(kLimit+1);
1267 <    eSin.resize(kLimit+1);
1266 >    elc.resize(kLimit+1);
1267 >    emc.resize(kLimit+1);
1268 >    enc.resize(kLimit+1);
1269 >    els.resize(kLimit+1);
1270 >    ems.resize(kLimit+1);
1271 >    ens.resize(kLimit+1);
1272 >
1273      for (int j = 0; j < kLimit+1; j++) {
1274 <      eCos[j].resize(nMax);
1275 <      eSin[j].resize(nMax);
1274 >      elc[j].resize(nMax);
1275 >      emc[j].resize(nMax);
1276 >      enc[j].resize(nMax);
1277 >      els[j].resize(nMax);
1278 >      ems[j].resize(nMax);
1279 >      ens[j].resize(nMax);
1280      }
1281      
1282      Vector3d t( 2.0 * M_PI );
1283      t.Vdiv(t, box);
1284  
1258 –    
1285      SimInfo::MoleculeIterator mi;
1286      Molecule::AtomIterator ai;
1287      int i;
1288      Vector3d r;
1289      Vector3d tt;
1264 –    Vector3d w;
1265 –    Vector3d u;
1266 –    Vector3d a;
1267 –    Vector3d b;
1290      
1291      for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1292           mol = info_->nextMolecule(mi)) {
# Line 1277 | Line 1299 | namespace OpenMD {
1299          
1300          tt.Vmul(t, r);
1301  
1302 <        
1303 <        eCos[1][i] = Vector3d(1.0, 1.0, 1.0);
1304 <        eSin[1][i] = Vector3d(0.0, 0.0, 0.0);
1305 <        eCos[2][i] = Vector3d(cos(tt.x()), cos(tt.y()), cos(tt.z()));
1306 <        eSin[2][i] = Vector3d(sin(tt.x()), sin(tt.y()), sin(tt.z()));
1302 >        elc[1][i] = 1.0;
1303 >        emc[1][i] = 1.0;
1304 >        enc[1][i] = 1.0;
1305 >        els[1][i] = 0.0;
1306 >        ems[1][i] = 0.0;
1307 >        ens[1][i] = 0.0;
1308  
1309 <        u = eCos[2][i];
1310 <        w = eSin[2][i];
1309 >        elc[2][i] = cos(tt.x());
1310 >        emc[2][i] = cos(tt.y());
1311 >        enc[2][i] = cos(tt.z());
1312 >        els[2][i] = sin(tt.x());
1313 >        ems[2][i] = sin(tt.y());
1314 >        ens[2][i] = sin(tt.z());
1315          
1316          for(int l = 3; l <= kLimit; l++) {
1317 <          eCos[l][i].x() = eCos[l-1][i].x()*eCos[2][i].x() - eSin[l-1][i].x()*eSin[2][i].x();
1318 <          eCos[l][i].y() = eCos[l-1][i].y()*eCos[2][i].y() - eSin[l-1][i].y()*eSin[2][i].y();
1319 <          eCos[l][i].z() = eCos[l-1][i].z()*eCos[2][i].z() - eSin[l-1][i].z()*eSin[2][i].z();
1320 <          
1321 <          eSin[l][i].x() = eSin[l-1][i].x()*eCos[2][i].x() + eCos[l-1][i].x()*eSin[2][i].x();
1322 <          eSin[l][i].y() = eSin[l-1][i].y()*eCos[2][i].y() + eCos[l-1][i].y()*eSin[2][i].y();
1296 <          eSin[l][i].z() = eSin[l-1][i].z()*eCos[2][i].z() + eCos[l-1][i].z()*eSin[2][i].z();
1297 <
1298 <
1299 <          // a.Vmul(eCos[l-1][i], u);
1300 <          // b.Vmul(eSin[l-1][i], w);
1301 <          // eCos[l][i] = a - b;
1302 <          // a.Vmul(eSin[l-1][i], u);
1303 <          // b.Vmul(eCos[l-1][i], w);
1304 <          // eSin[l][i] = a + b;
1305 <        
1317 >          elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i];
1318 >          emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i];
1319 >          enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i];
1320 >          els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i];
1321 >          ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i];
1322 >          ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i];
1323          }
1324        }
1325      }
# Line 1346 | Line 1363 | namespace OpenMD {
1363      std::vector<RealType> qks(nMax, 0.0);
1364      std::vector<Vector3d> dxk(nMax, V3Zero);
1365      std::vector<Vector3d> qxk(nMax, V3Zero);
1366 <    
1366 >    RealType rl, rm, rn;
1367 >    Vector3d kVec;
1368 >    Vector3d Qk;
1369 >    Mat3x3d k2;
1370 >    RealType ckcs, ckss, dkcs, dkss, qkcs, qkss;
1371 >    int atid;
1372 >    ElectrostaticAtomData data;
1373 >    RealType C, dk, qk;
1374 >    Vector3d D;
1375 >    Mat3x3d  Q;
1376 >
1377      int mMin = kLimit;
1378      int nMin = kLimit + 1;
1379      for (int l = 1; l <= kLimit; l++) {
1380        int ll = l - 1;
1381 <      RealType rl = xcl * float(ll);
1381 >      rl = xcl * float(ll);
1382        for (int mmm = mMin; mmm <= kLim2; mmm++) {
1383          int mm = mmm - kLimit;
1384          int m = abs(mm) + 1;
1385 <        RealType rm = ycl * float(mm);
1385 >        rm = ycl * float(mm);
1386          // Set temporary products of exponential terms
1387          for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1388               mol = info_->nextMolecule(mi)) {
# Line 1364 | Line 1391 | namespace OpenMD {
1391              
1392              i = atom->getLocalIndex();
1393              if(mm < 0) {
1394 <              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1395 <                     + eSin[l][i].x()*eSin[m][i].y();
1369 <              slm[i] = eSin[l][i].x()*eCos[m][i].y()
1370 <                     - eSin[m][i].y()*eCos[l][i].x();
1394 >              clm[i]=elc[l][i]*emc[m][i]+els[l][i]*ems[m][i];
1395 >              slm[i]=els[l][i]*emc[m][i]-ems[m][i]*elc[l][i];
1396              } else {
1397 <              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1398 <                     - eSin[l][i].x()*eSin[m][i].y();
1374 <              slm[i] = eSin[l][i].x()*eCos[m][i].y()
1375 <                     + eSin[m][i].y()*eCos[l][i].x();              
1397 >              clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i];
1398 >              slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i];
1399              }
1400            }
1401          }
1402          for (int nnn = nMin; nnn <= kLim2; nnn++) {
1403            int nn = nnn - kLimit;          
1404            int n = abs(nn) + 1;
1405 <          RealType rn = zcl * float(nn);
1405 >          rn = zcl * float(nn);
1406            // Test on magnitude of k vector:
1407            int kk=ll*ll + mm*mm + nn*nn;
1408            if(kk <= kSqLim) {
1409 <            Vector3d kVec = Vector3d(rl, rm, rn);
1410 <            Mat3x3d  k2 = outProduct(kVec, kVec);
1409 >            kVec = Vector3d(rl, rm, rn);
1410 >            k2 = outProduct(kVec, kVec);
1411              // Calculate exp(ikr) terms
1412              for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1413                   mol = info_->nextMolecule(mi)) {
# Line 1393 | Line 1416 | namespace OpenMD {
1416                  i = atom->getLocalIndex();
1417                  
1418                  if (nn < 0) {
1419 <                  ckr[i]=clm[i]*eCos[n][i].z()+slm[i]*eSin[n][i].z();
1420 <                  skr[i]=slm[i]*eCos[n][i].z()-clm[i]*eSin[n][i].z();
1419 >                  ckr[i]=clm[i]*enc[n][i]+slm[i]*ens[n][i];
1420 >                  skr[i]=slm[i]*enc[n][i]-clm[i]*ens[n][i];
1421 >
1422                  } else {
1423 <                  ckr[i]=clm[i]*eCos[n][i].z()-slm[i]*eSin[n][i].z();
1424 <                  skr[i]=slm[i]*eCos[n][i].z()+clm[i]*eSin[n][i].z();
1423 >                  ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i];
1424 >                  skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i];
1425                  }
1426                }
1427              }
# Line 1410 | Line 1434 | namespace OpenMD {
1434                    atom = mol->nextAtom(ai)) {
1435                  i = atom->getLocalIndex();
1436                  int atid = atom->getAtomType()->getIdent();
1437 <                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1437 >                data = ElectrostaticMap[Etids[atid]];
1438                                
1439                  if (data.is_Charge) {
1440 <                  RealType C = data.fixedCharge;
1440 >                  C = data.fixedCharge;
1441                    if (atom->isFluctuatingCharge()) C += atom->getFlucQPos();
1442                    ckc[i] = C * ckr[i];
1443                    cks[i] = C * skr[i];
1444                  }
1445                  
1446                  if (data.is_Dipole) {
1447 <                  Vector3d D = atom->getDipole() * mPoleConverter;
1448 <                  RealType dk = dot(D, kVec);
1447 >                  D = atom->getDipole() * mPoleConverter;
1448 >                  dk = dot(D, kVec);
1449                    dxk[i] = cross(D, kVec);
1450                    dkc[i] = dk * ckr[i];
1451                    dks[i] = dk * skr[i];
1452                  }
1453                  if (data.is_Quadrupole) {
1454 <                  Mat3x3d Q = atom->getQuadrupole();
1455 <                  Q *= mPoleConverter;
1456 <                  RealType qk = - doubleDot(Q, k2);
1457 <                  // RealType qk = -( Q * k2 ).trace();
1434 <                  qxk[i] = -2.0 * cross(k2, Q);
1454 >                  Q = atom->getQuadrupole() * mPoleConverter;
1455 >                  Qk = Q * kVec;                  
1456 >                  qk = dot(kVec, Qk);
1457 >                  qxk[i] = -cross(kVec, Qk);
1458                    qkc[i] = qk * ckr[i];
1459                    qks[i] = qk * skr[i];
1460                  }              
# Line 1440 | Line 1463 | namespace OpenMD {
1463  
1464              // calculate vector sums
1465              
1466 <            RealType ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1467 <            RealType ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1468 <            RealType dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1469 <            RealType dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1470 <            RealType qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1471 <            RealType qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1449 <
1466 >            ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1467 >            ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1468 >            dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1469 >            dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1470 >            qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1471 >            qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1472              
1473   #ifdef IS_MPI
1474 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE,
1475 <                                      MPI::SUM);
1476 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckss, 1, MPI::REALTYPE,
1477 <                                      MPI::SUM);
1478 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkcs, 1, MPI::REALTYPE,
1479 <                                      MPI::SUM);
1480 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkss, 1, MPI::REALTYPE,
1481 <                                      MPI::SUM);
1482 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkcs, 1, MPI::REALTYPE,
1483 <                                      MPI::SUM);
1484 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkss, 1, MPI::REALTYPE,
1485 <                                      MPI::SUM);
1474 >            MPI_Allreduce(MPI_IN_PLACE, &ckcs, 1, MPI_REALTYPE,
1475 >                          MPI_SUM, MPI_COMM_WORLD);
1476 >            MPI_Allreduce(MPI_IN_PLACE, &ckss, 1, MPI_REALTYPE,
1477 >                          MPI_SUM, MPI_COMM_WORLD);
1478 >            MPI_Allreduce(MPI_IN_PLACE, &dkcs, 1, MPI_REALTYPE,
1479 >                          MPI_SUM, MPI_COMM_WORLD);
1480 >            MPI_Allreduce(MPI_IN_PLACE, &dkss, 1, MPI_REALTYPE,
1481 >                          MPI_SUM, MPI_COMM_WORLD);
1482 >            MPI_Allreduce(MPI_IN_PLACE, &qkcs, 1, MPI_REALTYPE,
1483 >                          MPI_SUM, MPI_COMM_WORLD);
1484 >            MPI_Allreduce(MPI_IN_PLACE, &qkss, 1, MPI_REALTYPE,
1485 >                          MPI_SUM, MPI_COMM_WORLD);
1486   #endif        
1487              
1488              // Accumulate potential energy and virial contribution:
1489  
1490 <            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1491 <                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkss));
1490 >            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1491 >                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs));
1492  
1493 <            kVir -= 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss
1494 <                                         +4.0*(ckss*dkcs-ckcs*dkss)
1495 <                                         +3.0*(dkcs*dkcs+dkss*dkss)
1496 <                                         -6.0*(ckss*qkss+ckcs*qkcs)
1497 <                                         +8.0*(dkss*qkcs-dkcs*qkss)
1498 <                                         +5.0*(qkss*qkss+qkcs*qkcs));
1493 >            kVir += 2.0 * rvol  * AK[kk]*(ckcs*ckcs+ckss*ckss
1494 >                                          +4.0*(ckss*dkcs-ckcs*dkss)
1495 >                                          +3.0*(dkcs*dkcs+dkss*dkss)
1496 >                                          -6.0*(ckss*qkss+ckcs*qkcs)
1497 >                                          +8.0*(dkss*qkcs-dkcs*qkss)
1498 >                                          +5.0*(qkss*qkss+qkcs*qkcs));
1499              
1500              // Calculate force and torque for each site:
1501              
# Line 1483 | Line 1505 | namespace OpenMD {
1505                    atom = mol->nextAtom(ai)) {
1506                  
1507                  i = atom->getLocalIndex();
1508 <                int atid = atom->getAtomType()->getIdent();
1509 <                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1510 <                
1508 >                atid = atom->getAtomType()->getIdent();
1509 >                data = ElectrostaticMap[Etids[atid]];
1510 >
1511                  RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs)
1512                                       - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss));
1513                  RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
1514                                           -ckr[i]*(ckss+dkcs-qkss));
1515 <                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)+
1516 <                                             skr[i]*(ckss+dkcs-qkss));
1515 >                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
1516 >                                            +skr[i]*(ckss+dkcs-qkss));
1517                
1518                  atom->addFrc( 4.0 * rvol * qfrc * kVec );
1519 <                
1519 >
1520 >                if (atom->isFluctuatingCharge()) {
1521 >                  atom->addFlucQFrc( - 2.0 * rvol * qtrq2 );
1522 >                }
1523 >                  
1524                  if (data.is_Dipole) {
1525                    atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] );
1526                  }
# Line 1509 | Line 1535 | namespace OpenMD {
1535        }
1536        mMin = 1;
1537      }
1538 <    cerr << "kPot = " << kPot << "\n";
1513 <    pot[ELECTROSTATIC_FAMILY] += kPot;  
1538 >    pot += kPot;  
1539    }
1540   }

Diff Legend

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