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 1994 by gezelter, Wed Apr 30 18:50:45 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 752 | 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 764 | Line 782 | namespace OpenMD {
782      // Excluded potential that is still computed for fluctuating charges
783      excluded_Pot= 0.0;
784  
767 –
785      // some variables we'll need independent of electrostatic type:
786  
787      ri = 1.0 /  *(idat.rij);
# Line 827 | 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 836 | 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 847 | 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 863 | 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 870 | 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 881 | 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 <    
917 >        
918 >
919      if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
920        J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
921      }    
922 <    
922 >
923      if (a_is_Charge) {    
924        
925        if (b_is_Charge) {
926          pref =  pre11_ * *(idat.electroMult);      
927          U  += C_a * C_b * pref * v01;
928          F  += C_a * C_b * pref * dv01 * rhat;
929 <        
929 >
930          // If this is an excluded pair, there are still indirect
931          // interactions via the reaction field we must worry about:
932  
# Line 905 | Line 935 | namespace OpenMD {
935            indirect_Pot += rfContrib;
936            indirect_F   += rfContrib * 2.0 * ri * rhat;
937          }
938 <        
938 >
939          // Fluctuating charge forces are handled via Coulomb integrals
940          // for excluded pairs (i.e. those connected via bonds) and
941          // with the standard charge-charge interaction otherwise.
942  
943 <        if (idat.excluded) {          
943 >        if (idat.excluded) {
944            if (a_is_Fluctuating || b_is_Fluctuating) {
945              coulInt = J->getValueAt( *(idat.rij) );
946 <            if (a_is_Fluctuating)  dUdCa += coulInt * C_b;
947 <            if (b_is_Fluctuating)  dUdCb += coulInt * C_a;
948 <            excluded_Pot += C_a * C_b * coulInt;
919 <          }          
946 >            if (a_is_Fluctuating) dUdCa += C_b * coulInt;
947 >            if (b_is_Fluctuating) dUdCb += C_a * coulInt;          
948 >          }
949          } else {
950            if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
951 <          if (a_is_Fluctuating) dUdCb += C_a * pref * v01;
952 <        }
951 >          if (b_is_Fluctuating) dUdCb += C_a * pref * v01;
952 >        }              
953        }
954  
955        if (b_is_Dipole) {
# Line 986 | Line 1015 | namespace OpenMD {
1015          F  -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat;
1016          Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
1017          Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
989 –
1018          // Even if we excluded this pair from direct interactions, we
1019          // still have the reaction-field-mediated dipole-dipole
1020          // interaction:
# Line 1046 | Line 1074 | namespace OpenMD {
1074          trQaQb = QaQb.trace();
1075          rQaQb = rhat * QaQb;
1076          QaQbr = QaQb * rhat;
1077 <        QaxQb = cross(Q_a, Q_b);
1077 >        QaxQb = mCross(Q_a, Q_b);
1078          rQaQbr = dot(rQa, Qbr);
1079          rQaxQbr = cross(rQa, Qbr);
1080          
# Line 1077 | Line 1105 | namespace OpenMD {
1105          //             + 4.0 * cross(rhat, QbQar)
1106  
1107          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1080 –
1108        }
1109      }
1110  
# Line 1086 | Line 1113 | namespace OpenMD {
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);
1122      if (b_is_Fluctuating) *(idat.dVdFQ2) += dUdCb * *(idat.sw);
1123  
# Line 1140 | Line 1172 | namespace OpenMD {
1172          
1173      if (i_is_Fluctuating) {
1174        C_a += *(sdat.flucQ);
1175 <      // dVdFQ is really a force, so this is negative the derivative
1176 <      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1177 <      (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1178 <        (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1175 >
1176 >      flucQ_->getSelfInteraction(sdat.atid, *(sdat.flucQ),  
1177 >                                 (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY],
1178 >                                 *(sdat.flucQfrc) );
1179 >
1180      }
1181  
1182      switch (summationMethod_) {
# Line 1194 | Line 1227 | namespace OpenMD {
1227    }
1228  
1229  
1230 <  void Electrostatic::ReciprocalSpaceSum(potVec& pot) {
1230 >  void Electrostatic::ReciprocalSpaceSum(RealType& pot) {
1231      
1232      RealType kPot = 0.0;
1233      RealType kVir = 0.0;
# Line 1240 | Line 1273 | namespace OpenMD {
1273      
1274      // Calculate and store exponential factors  
1275      
1276 <    vector<vector<Vector3d> > eCos;
1277 <    vector<vector<Vector3d> > eSin;
1276 >    vector<vector<RealType> > elc;
1277 >    vector<vector<RealType> > emc;
1278 >    vector<vector<RealType> > enc;
1279 >    vector<vector<RealType> > els;
1280 >    vector<vector<RealType> > ems;
1281 >    vector<vector<RealType> > ens;
1282      
1283      int nMax = info_->getNAtoms();
1284      
1285 <    eCos.resize(kLimit+1);
1286 <    eSin.resize(kLimit+1);
1285 >    elc.resize(kLimit+1);
1286 >    emc.resize(kLimit+1);
1287 >    enc.resize(kLimit+1);
1288 >    els.resize(kLimit+1);
1289 >    ems.resize(kLimit+1);
1290 >    ens.resize(kLimit+1);
1291 >
1292      for (int j = 0; j < kLimit+1; j++) {
1293 <      eCos[j].resize(nMax);
1294 <      eSin[j].resize(nMax);
1293 >      elc[j].resize(nMax);
1294 >      emc[j].resize(nMax);
1295 >      enc[j].resize(nMax);
1296 >      els[j].resize(nMax);
1297 >      ems[j].resize(nMax);
1298 >      ens[j].resize(nMax);
1299      }
1300      
1301      Vector3d t( 2.0 * M_PI );
1302      t.Vdiv(t, box);
1303  
1258 –    
1304      SimInfo::MoleculeIterator mi;
1305      Molecule::AtomIterator ai;
1306      int i;
1307      Vector3d r;
1308      Vector3d tt;
1264 –    Vector3d w;
1265 –    Vector3d u;
1266 –    Vector3d a;
1267 –    Vector3d b;
1309      
1310      for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1311           mol = info_->nextMolecule(mi)) {
# Line 1277 | Line 1318 | namespace OpenMD {
1318          
1319          tt.Vmul(t, r);
1320  
1321 <        
1322 <        eCos[1][i] = Vector3d(1.0, 1.0, 1.0);
1323 <        eSin[1][i] = Vector3d(0.0, 0.0, 0.0);
1324 <        eCos[2][i] = Vector3d(cos(tt.x()), cos(tt.y()), cos(tt.z()));
1325 <        eSin[2][i] = Vector3d(sin(tt.x()), sin(tt.y()), sin(tt.z()));
1321 >        elc[1][i] = 1.0;
1322 >        emc[1][i] = 1.0;
1323 >        enc[1][i] = 1.0;
1324 >        els[1][i] = 0.0;
1325 >        ems[1][i] = 0.0;
1326 >        ens[1][i] = 0.0;
1327  
1328 <        u = eCos[2][i];
1329 <        w = eSin[2][i];
1328 >        elc[2][i] = cos(tt.x());
1329 >        emc[2][i] = cos(tt.y());
1330 >        enc[2][i] = cos(tt.z());
1331 >        els[2][i] = sin(tt.x());
1332 >        ems[2][i] = sin(tt.y());
1333 >        ens[2][i] = sin(tt.z());
1334          
1335          for(int l = 3; l <= kLimit; l++) {
1336 <          eCos[l][i].x() = eCos[l-1][i].x()*eCos[2][i].x() - eSin[l-1][i].x()*eSin[2][i].x();
1337 <          eCos[l][i].y() = eCos[l-1][i].y()*eCos[2][i].y() - eSin[l-1][i].y()*eSin[2][i].y();
1338 <          eCos[l][i].z() = eCos[l-1][i].z()*eCos[2][i].z() - eSin[l-1][i].z()*eSin[2][i].z();
1339 <          
1340 <          eSin[l][i].x() = eSin[l-1][i].x()*eCos[2][i].x() + eCos[l-1][i].x()*eSin[2][i].x();
1341 <          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 <        
1336 >          elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i];
1337 >          emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i];
1338 >          enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i];
1339 >          els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i];
1340 >          ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i];
1341 >          ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i];
1342          }
1343        }
1344      }
# Line 1346 | Line 1382 | namespace OpenMD {
1382      std::vector<RealType> qks(nMax, 0.0);
1383      std::vector<Vector3d> dxk(nMax, V3Zero);
1384      std::vector<Vector3d> qxk(nMax, V3Zero);
1385 <    
1385 >    RealType rl, rm, rn;
1386 >    Vector3d kVec;
1387 >    Vector3d Qk;
1388 >    Mat3x3d k2;
1389 >    RealType ckcs, ckss, dkcs, dkss, qkcs, qkss;
1390 >    int atid;
1391 >    ElectrostaticAtomData data;
1392 >    RealType C, dk, qk;
1393 >    Vector3d D;
1394 >    Mat3x3d  Q;
1395 >
1396      int mMin = kLimit;
1397      int nMin = kLimit + 1;
1398      for (int l = 1; l <= kLimit; l++) {
1399        int ll = l - 1;
1400 <      RealType rl = xcl * float(ll);
1400 >      rl = xcl * float(ll);
1401        for (int mmm = mMin; mmm <= kLim2; mmm++) {
1402          int mm = mmm - kLimit;
1403          int m = abs(mm) + 1;
1404 <        RealType rm = ycl * float(mm);
1404 >        rm = ycl * float(mm);
1405          // Set temporary products of exponential terms
1406          for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1407               mol = info_->nextMolecule(mi)) {
# Line 1364 | Line 1410 | namespace OpenMD {
1410              
1411              i = atom->getLocalIndex();
1412              if(mm < 0) {
1413 <              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1414 <                     + 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();
1413 >              clm[i]=elc[l][i]*emc[m][i]+els[l][i]*ems[m][i];
1414 >              slm[i]=els[l][i]*emc[m][i]-ems[m][i]*elc[l][i];
1415              } else {
1416 <              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1417 <                     - 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();              
1416 >              clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i];
1417 >              slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i];
1418              }
1419            }
1420          }
1421          for (int nnn = nMin; nnn <= kLim2; nnn++) {
1422            int nn = nnn - kLimit;          
1423            int n = abs(nn) + 1;
1424 <          RealType rn = zcl * float(nn);
1424 >          rn = zcl * float(nn);
1425            // Test on magnitude of k vector:
1426            int kk=ll*ll + mm*mm + nn*nn;
1427            if(kk <= kSqLim) {
1428 <            Vector3d kVec = Vector3d(rl, rm, rn);
1429 <            Mat3x3d  k2 = outProduct(kVec, kVec);
1428 >            kVec = Vector3d(rl, rm, rn);
1429 >            k2 = outProduct(kVec, kVec);
1430              // Calculate exp(ikr) terms
1431              for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1432                   mol = info_->nextMolecule(mi)) {
# Line 1393 | Line 1435 | namespace OpenMD {
1435                  i = atom->getLocalIndex();
1436                  
1437                  if (nn < 0) {
1438 <                  ckr[i]=clm[i]*eCos[n][i].z()+slm[i]*eSin[n][i].z();
1439 <                  skr[i]=slm[i]*eCos[n][i].z()-clm[i]*eSin[n][i].z();
1438 >                  ckr[i]=clm[i]*enc[n][i]+slm[i]*ens[n][i];
1439 >                  skr[i]=slm[i]*enc[n][i]-clm[i]*ens[n][i];
1440 >
1441                  } else {
1442 <                  ckr[i]=clm[i]*eCos[n][i].z()-slm[i]*eSin[n][i].z();
1443 <                  skr[i]=slm[i]*eCos[n][i].z()+clm[i]*eSin[n][i].z();
1442 >                  ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i];
1443 >                  skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i];
1444                  }
1445                }
1446              }
# Line 1410 | Line 1453 | namespace OpenMD {
1453                    atom = mol->nextAtom(ai)) {
1454                  i = atom->getLocalIndex();
1455                  int atid = atom->getAtomType()->getIdent();
1456 <                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1456 >                data = ElectrostaticMap[Etids[atid]];
1457                                
1458                  if (data.is_Charge) {
1459 <                  RealType C = data.fixedCharge;
1459 >                  C = data.fixedCharge;
1460                    if (atom->isFluctuatingCharge()) C += atom->getFlucQPos();
1461                    ckc[i] = C * ckr[i];
1462                    cks[i] = C * skr[i];
1463                  }
1464                  
1465                  if (data.is_Dipole) {
1466 <                  Vector3d D = atom->getDipole() * mPoleConverter;
1467 <                  RealType dk = dot(D, kVec);
1466 >                  D = atom->getDipole() * mPoleConverter;
1467 >                  dk = dot(D, kVec);
1468                    dxk[i] = cross(D, kVec);
1469                    dkc[i] = dk * ckr[i];
1470                    dks[i] = dk * skr[i];
1471                  }
1472                  if (data.is_Quadrupole) {
1473 <                  Mat3x3d Q = atom->getQuadrupole();
1474 <                  Q *= mPoleConverter;
1475 <                  RealType qk = - doubleDot(Q, k2);
1476 <                  // RealType qk = -( Q * k2 ).trace();
1434 <                  qxk[i] = -2.0 * cross(k2, Q);
1473 >                  Q = atom->getQuadrupole() * mPoleConverter;
1474 >                  Qk = Q * kVec;                  
1475 >                  qk = dot(kVec, Qk);
1476 >                  qxk[i] = -cross(kVec, Qk);
1477                    qkc[i] = qk * ckr[i];
1478                    qks[i] = qk * skr[i];
1479                  }              
# Line 1440 | Line 1482 | namespace OpenMD {
1482  
1483              // calculate vector sums
1484              
1485 <            RealType ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1486 <            RealType ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1487 <            RealType dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1488 <            RealType dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1489 <            RealType qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1490 <            RealType qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1449 <
1485 >            ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1486 >            ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1487 >            dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1488 >            dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1489 >            qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1490 >            qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1491              
1492   #ifdef IS_MPI
1493 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE,
1494 <                                      MPI::SUM);
1495 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckss, 1, MPI::REALTYPE,
1496 <                                      MPI::SUM);
1497 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkcs, 1, MPI::REALTYPE,
1498 <                                      MPI::SUM);
1499 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkss, 1, MPI::REALTYPE,
1500 <                                      MPI::SUM);
1501 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkcs, 1, MPI::REALTYPE,
1502 <                                      MPI::SUM);
1503 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkss, 1, MPI::REALTYPE,
1504 <                                      MPI::SUM);
1493 >            MPI_Allreduce(MPI_IN_PLACE, &ckcs, 1, MPI_REALTYPE,
1494 >                          MPI_SUM, MPI_COMM_WORLD);
1495 >            MPI_Allreduce(MPI_IN_PLACE, &ckss, 1, MPI_REALTYPE,
1496 >                          MPI_SUM, MPI_COMM_WORLD);
1497 >            MPI_Allreduce(MPI_IN_PLACE, &dkcs, 1, MPI_REALTYPE,
1498 >                          MPI_SUM, MPI_COMM_WORLD);
1499 >            MPI_Allreduce(MPI_IN_PLACE, &dkss, 1, MPI_REALTYPE,
1500 >                          MPI_SUM, MPI_COMM_WORLD);
1501 >            MPI_Allreduce(MPI_IN_PLACE, &qkcs, 1, MPI_REALTYPE,
1502 >                          MPI_SUM, MPI_COMM_WORLD);
1503 >            MPI_Allreduce(MPI_IN_PLACE, &qkss, 1, MPI_REALTYPE,
1504 >                          MPI_SUM, MPI_COMM_WORLD);
1505   #endif        
1506              
1507              // Accumulate potential energy and virial contribution:
1508  
1509 <            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1510 <                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkss));
1509 >            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1510 >                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs));
1511  
1512 <            kVir -= 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss
1513 <                                         +4.0*(ckss*dkcs-ckcs*dkss)
1514 <                                         +3.0*(dkcs*dkcs+dkss*dkss)
1515 <                                         -6.0*(ckss*qkss+ckcs*qkcs)
1516 <                                         +8.0*(dkss*qkcs-dkcs*qkss)
1517 <                                         +5.0*(qkss*qkss+qkcs*qkcs));
1512 >            kVir += 2.0 * rvol  * AK[kk]*(ckcs*ckcs+ckss*ckss
1513 >                                          +4.0*(ckss*dkcs-ckcs*dkss)
1514 >                                          +3.0*(dkcs*dkcs+dkss*dkss)
1515 >                                          -6.0*(ckss*qkss+ckcs*qkcs)
1516 >                                          +8.0*(dkss*qkcs-dkcs*qkss)
1517 >                                          +5.0*(qkss*qkss+qkcs*qkcs));
1518              
1519              // Calculate force and torque for each site:
1520              
# Line 1483 | Line 1524 | namespace OpenMD {
1524                    atom = mol->nextAtom(ai)) {
1525                  
1526                  i = atom->getLocalIndex();
1527 <                int atid = atom->getAtomType()->getIdent();
1528 <                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1529 <                
1527 >                atid = atom->getAtomType()->getIdent();
1528 >                data = ElectrostaticMap[Etids[atid]];
1529 >
1530                  RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs)
1531                                       - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss));
1532                  RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
1533                                           -ckr[i]*(ckss+dkcs-qkss));
1534 <                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)+
1535 <                                             skr[i]*(ckss+dkcs-qkss));
1534 >                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
1535 >                                            +skr[i]*(ckss+dkcs-qkss));
1536                
1537                  atom->addFrc( 4.0 * rvol * qfrc * kVec );
1538 <                
1538 >
1539 >                if (atom->isFluctuatingCharge()) {
1540 >                  atom->addFlucQFrc( - 2.0 * rvol * qtrq2 );
1541 >                }
1542 >                  
1543                  if (data.is_Dipole) {
1544                    atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] );
1545                  }
# Line 1509 | Line 1554 | namespace OpenMD {
1554        }
1555        mMin = 1;
1556      }
1557 <    cerr << "kPot = " << kPot << "\n";
1513 <    pot[ELECTROSTATIC_FAMILY] += kPot;  
1557 >    pot += kPot;  
1558    }
1559 +
1560 +  void Electrostatic::getSitePotentials(Atom* a1, Atom* a2, bool excluded,
1561 +                                        RealType &spot1, RealType &spot2) {
1562 +
1563 +    if (!initialized_) {
1564 +      cerr << "initializing\n";
1565 +      initialize();
1566 +      cerr << "done\n";
1567 +    }
1568 +
1569 +    const RealType mPoleConverter = 0.20819434;
1570 +
1571 +    AtomType* atype1 = a1->getAtomType();
1572 +    AtomType* atype2 = a2->getAtomType();
1573 +    int atid1 = atype1->getIdent();
1574 +    int atid2 = atype2->getIdent();
1575 +    data1 = ElectrostaticMap[Etids[atid1]];
1576 +    data2 = ElectrostaticMap[Etids[atid2]];
1577 +
1578 +    Pa = 0.0;  // Site potential at site a
1579 +    Pb = 0.0;  // Site potential at site b
1580 +
1581 +    Vector3d d = a2->getPos() - a1->getPos();
1582 +    info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(d);
1583 +    RealType rij = d.length();
1584 +    // some variables we'll need independent of electrostatic type:
1585 +
1586 +    RealType ri = 1.0 /  rij;
1587 +    rhat = d  * ri;
1588 +      
1589 +
1590 +    if ((rij >= cutoffRadius_) || excluded) {
1591 +      spot1 = 0.0;
1592 +      spot2 = 0.0;
1593 +      return;
1594 +    }
1595 +
1596 +    // logicals
1597 +
1598 +    a_is_Charge = data1.is_Charge;
1599 +    a_is_Dipole = data1.is_Dipole;
1600 +    a_is_Quadrupole = data1.is_Quadrupole;
1601 +    a_is_Fluctuating = data1.is_Fluctuating;
1602 +
1603 +    b_is_Charge = data2.is_Charge;
1604 +    b_is_Dipole = data2.is_Dipole;
1605 +    b_is_Quadrupole = data2.is_Quadrupole;
1606 +    b_is_Fluctuating = data2.is_Fluctuating;
1607 +
1608 +    // Obtain all of the required radial function values from the
1609 +    // spline structures:
1610 +    
1611 +
1612 +    if (a_is_Charge || b_is_Charge) {
1613 +      v01 = v01s->getValueAt(rij);
1614 +    }
1615 +    if (a_is_Dipole || b_is_Dipole) {
1616 +      v11 = v11s->getValueAt(rij);
1617 +      v11or = ri * v11;
1618 +    }
1619 +    if (a_is_Quadrupole || b_is_Quadrupole) {
1620 +      v21 = v21s->getValueAt(rij);
1621 +      v22 = v22s->getValueAt(rij);
1622 +      v22or = ri * v22;
1623 +    }      
1624 +
1625 +    if (a_is_Charge) {
1626 +      C_a = data1.fixedCharge;
1627 +      
1628 +      if (a_is_Fluctuating) {
1629 +        C_a += a1->getFlucQPos();
1630 +      }
1631 +      
1632 +      Pb += C_a *  pre11_ * v01;      
1633 +    }
1634 +    
1635 +    if (a_is_Dipole) {
1636 +      D_a = a1->getDipole() * mPoleConverter;
1637 +      rdDa = dot(rhat, D_a);
1638 +      Pb +=  pre12_ * v11 * rdDa;      
1639 +    }
1640 +    
1641 +    if (a_is_Quadrupole) {
1642 +      Q_a = a1->getQuadrupole() * mPoleConverter;
1643 +      trQa =  Q_a.trace();
1644 +      Qar =   Q_a * rhat;
1645 +      rdQar = dot(rhat, Qar);
1646 +      Pb += pre14_ * (v21 * trQa + v22 * rdQar);      
1647 +    }
1648 +    
1649 +    if (b_is_Charge) {
1650 +      C_b = data2.fixedCharge;
1651 +      
1652 +      if (b_is_Fluctuating)
1653 +        C_b += a2->getFlucQPos();
1654 +      
1655 +      Pa += C_b *  pre11_ * v01;
1656 +    }
1657 +    
1658 +    if (b_is_Dipole) {
1659 +      D_b = a2->getDipole() * mPoleConverter;
1660 +      rdDb = dot(rhat, D_b);
1661 +      Pa += pre12_ * v11 * rdDb;
1662 +    }
1663 +    
1664 +    if (b_is_Quadrupole) {
1665 +      Q_a = a2->getQuadrupole() * mPoleConverter;
1666 +      trQb =  Q_b.trace();
1667 +      Qbr =   Q_b * rhat;
1668 +      rdQbr = dot(rhat, Qbr);
1669 +      Pa += pre14_ * (v21 * trQb + v22 * rdQbr);      
1670 +    }
1671 +
1672 +    spot1 = Pa;
1673 +    spot2 = Pb;
1674 +  }
1675   }

Diff Legend

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