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 2071 by gezelter, Sat Mar 7 21:41:51 2015 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),
65 –                                  forceField_(NULL), info_(NULL),
69                                    haveCutoffRadius_(false),
70                                    haveDampingAlpha_(false),
71                                    haveDielectric_(false),
72 <                                  haveElectroSplines_(false)
73 <  {}
72 >                                  haveElectroSplines_(false),
73 >                                  info_(NULL), forceField_(NULL)
74 >                                  
75 >  {
76 >    flucQ_ = new FluctuatingChargeForces(info_);
77 >  }
78    
79 +  void Electrostatic::setForceField(ForceField *ff) {
80 +    forceField_ = ff;
81 +    flucQ_->setForceField(forceField_);
82 +  }
83 +
84 +  void Electrostatic::setSimulatedAtomTypes(set<AtomType*> &simtypes) {
85 +    simTypes_ = simtypes;
86 +    flucQ_->setSimulatedAtomTypes(simTypes_);
87 +  }
88 +
89    void Electrostatic::initialize() {
90      
91      Globals* simParams_ = info_->getSimParams();
# Line 247 | Line 264 | namespace OpenMD {
264      
265      RealType b0c, b1c, b2c, b3c, b4c, b5c;
266      RealType db0c_1, db0c_2, db0c_3, db0c_4, db0c_5;
267 <    RealType a2, expTerm, invArootPi;
267 >    RealType a2, expTerm, invArootPi(0.0);
268      
269      RealType r = cutoffRadius_;
270      RealType r2 = r * r;
# Line 752 | Line 769 | namespace OpenMD {
769      Tb.zero(); // Torque on site b
770      Ea.zero(); // Electric field at site a
771      Eb.zero(); // Electric field at site b
772 +    Pa = 0.0;  // Site potential at site a
773 +    Pb = 0.0;  // Site potential at site b
774      dUdCa = 0.0; // fluctuating charge force at site a
775      dUdCb = 0.0; // fluctuating charge force at site a
776      
# Line 764 | Line 783 | namespace OpenMD {
783      // Excluded potential that is still computed for fluctuating charges
784      excluded_Pot= 0.0;
785  
767 –
786      // some variables we'll need independent of electrostatic type:
787  
788      ri = 1.0 /  *(idat.rij);
# Line 827 | Line 845 | namespace OpenMD {
845        if (idat.excluded) {
846          *(idat.skippedCharge2) += C_a;
847        } else {
848 <        // only do the field if we're not excluded:
848 >        // only do the field and site potentials if we're not excluded:
849          Eb -= C_a *  pre11_ * dv01 * rhat;
850 +        Pb += C_a *  pre11_ * v01;
851        }
852      }
853      
# Line 836 | Line 855 | namespace OpenMD {
855        D_a = *(idat.dipole1);
856        rdDa = dot(rhat, D_a);
857        rxDa = cross(rhat, D_a);
858 <      if (!idat.excluded)
858 >      if (!idat.excluded) {
859          Eb -=  pre12_ * ((dv11-v11or) * rdDa * rhat + v11or * D_a);
860 +        Pb +=  pre12_ * v11 * rdDa;
861 +      }
862 +
863      }
864      
865      if (a_is_Quadrupole) {
# Line 847 | Line 869 | namespace OpenMD {
869        rQa = rhat * Q_a;
870        rdQar = dot(rhat, Qar);
871        rxQar = cross(rhat, Qar);
872 <      if (!idat.excluded)
872 >      if (!idat.excluded) {
873          Eb -= pre14_ * (trQa * rhat * dv21 + 2.0 * Qar * v22or
874                          + rdQar * rhat * (dv22 - 2.0*v22or));
875 +        Pb += pre14_ * (v21 * trQa + v22 * rdQar);
876 +      }
877      }
878      
879      if (b_is_Charge) {
# Line 863 | Line 887 | namespace OpenMD {
887        } else {
888          // only do the field if we're not excluded:
889          Ea += C_b *  pre11_ * dv01 * rhat;
890 +        Pa += C_b *  pre11_ * v01;
891 +
892        }
893      }
894      
# Line 870 | Line 896 | namespace OpenMD {
896        D_b = *(idat.dipole2);
897        rdDb = dot(rhat, D_b);
898        rxDb = cross(rhat, D_b);
899 <      if (!idat.excluded)
899 >      if (!idat.excluded) {
900          Ea += pre12_ * ((dv11-v11or) * rdDb * rhat + v11or * D_b);
901 +        Pa += pre12_ * v11 * rdDb;
902 +      }
903      }
904      
905      if (b_is_Quadrupole) {
# Line 881 | Line 909 | namespace OpenMD {
909        rQb = rhat * Q_b;
910        rdQbr = dot(rhat, Qbr);
911        rxQbr = cross(rhat, Qbr);
912 <      if (!idat.excluded)
912 >      if (!idat.excluded) {
913          Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or
914                          + rdQbr * rhat * (dv22 - 2.0*v22or));
915 +        Pa += pre14_ * (v21 * trQb + v22 * rdQbr);
916 +      }
917      }
918 <    
918 >        
919 >
920      if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
921        J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
922      }    
923 <    
923 >
924      if (a_is_Charge) {    
925        
926        if (b_is_Charge) {
927          pref =  pre11_ * *(idat.electroMult);      
928          U  += C_a * C_b * pref * v01;
929          F  += C_a * C_b * pref * dv01 * rhat;
930 <        
930 >
931          // If this is an excluded pair, there are still indirect
932          // interactions via the reaction field we must worry about:
933  
# Line 905 | Line 936 | namespace OpenMD {
936            indirect_Pot += rfContrib;
937            indirect_F   += rfContrib * 2.0 * ri * rhat;
938          }
939 <        
939 >
940          // Fluctuating charge forces are handled via Coulomb integrals
941          // for excluded pairs (i.e. those connected via bonds) and
942          // with the standard charge-charge interaction otherwise.
943  
944 <        if (idat.excluded) {          
944 >        if (idat.excluded) {
945            if (a_is_Fluctuating || b_is_Fluctuating) {
946              coulInt = J->getValueAt( *(idat.rij) );
947 <            if (a_is_Fluctuating)  dUdCa += coulInt * C_b;
948 <            if (b_is_Fluctuating)  dUdCb += coulInt * C_a;
949 <            excluded_Pot += C_a * C_b * coulInt;
919 <          }          
947 >            if (a_is_Fluctuating) dUdCa += C_b * coulInt;
948 >            if (b_is_Fluctuating) dUdCb += C_a * coulInt;          
949 >          }
950          } else {
951            if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
952 <          if (a_is_Fluctuating) dUdCb += C_a * pref * v01;
953 <        }
952 >          if (b_is_Fluctuating) dUdCb += C_a * pref * v01;
953 >        }              
954        }
955  
956        if (b_is_Dipole) {
# Line 986 | Line 1016 | namespace OpenMD {
1016          F  -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat;
1017          Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
1018          Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
989 –
1019          // Even if we excluded this pair from direct interactions, we
1020          // still have the reaction-field-mediated dipole-dipole
1021          // interaction:
# Line 1046 | Line 1075 | namespace OpenMD {
1075          trQaQb = QaQb.trace();
1076          rQaQb = rhat * QaQb;
1077          QaQbr = QaQb * rhat;
1078 <        QaxQb = cross(Q_a, Q_b);
1078 >        QaxQb = mCross(Q_a, Q_b);
1079          rQaQbr = dot(rQa, Qbr);
1080          rQaxQbr = cross(rQa, Qbr);
1081          
# Line 1077 | Line 1106 | namespace OpenMD {
1106          //             + 4.0 * cross(rhat, QbQar)
1107  
1108          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1080 –
1109        }
1110      }
1111  
1112      if (idat.doElectricField) {
1113        *(idat.eField1) += Ea * *(idat.electroMult);
1114        *(idat.eField2) += Eb * *(idat.electroMult);
1115 +    }
1116 +
1117 +    if (idat.doSitePotential) {
1118 +      *(idat.sPot1) += Pa * *(idat.electroMult);
1119 +      *(idat.sPot2) += Pb * *(idat.electroMult);
1120      }
1121  
1122      if (a_is_Fluctuating) *(idat.dVdFQ1) += dUdCa * *(idat.sw);
# Line 1132 | Line 1165 | namespace OpenMD {
1165      bool i_is_Quadrupole = data.is_Quadrupole;
1166      bool i_is_Fluctuating = data.is_Fluctuating;
1167      RealType C_a = data.fixedCharge;  
1168 <    RealType self(0.0), preVal, DdD, trQ, trQQ;
1168 >    RealType self(0.0), preVal, DdD(0.0), trQ, trQQ;
1169  
1170      if (i_is_Dipole) {
1171        DdD = data.dipole.lengthSquare();
# Line 1140 | Line 1173 | namespace OpenMD {
1173          
1174      if (i_is_Fluctuating) {
1175        C_a += *(sdat.flucQ);
1176 <      // dVdFQ is really a force, so this is negative the derivative
1177 <      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1178 <      (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1179 <        (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1176 >
1177 >      flucQ_->getSelfInteraction(sdat.atid, *(sdat.flucQ),  
1178 >                                 (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY],
1179 >                                 *(sdat.flucQfrc) );
1180 >
1181      }
1182  
1183      switch (summationMethod_) {
# Line 1194 | Line 1228 | namespace OpenMD {
1228    }
1229  
1230  
1231 <  void Electrostatic::ReciprocalSpaceSum(potVec& pot) {
1231 >  void Electrostatic::ReciprocalSpaceSum(RealType& pot) {
1232      
1233      RealType kPot = 0.0;
1234      RealType kVir = 0.0;
# Line 1240 | Line 1274 | namespace OpenMD {
1274      
1275      // Calculate and store exponential factors  
1276      
1277 <    vector<vector<Vector3d> > eCos;
1278 <    vector<vector<Vector3d> > eSin;
1277 >    vector<vector<RealType> > elc;
1278 >    vector<vector<RealType> > emc;
1279 >    vector<vector<RealType> > enc;
1280 >    vector<vector<RealType> > els;
1281 >    vector<vector<RealType> > ems;
1282 >    vector<vector<RealType> > ens;
1283      
1284      int nMax = info_->getNAtoms();
1285      
1286 <    eCos.resize(kLimit+1);
1287 <    eSin.resize(kLimit+1);
1286 >    elc.resize(kLimit+1);
1287 >    emc.resize(kLimit+1);
1288 >    enc.resize(kLimit+1);
1289 >    els.resize(kLimit+1);
1290 >    ems.resize(kLimit+1);
1291 >    ens.resize(kLimit+1);
1292 >
1293      for (int j = 0; j < kLimit+1; j++) {
1294 <      eCos[j].resize(nMax);
1295 <      eSin[j].resize(nMax);
1294 >      elc[j].resize(nMax);
1295 >      emc[j].resize(nMax);
1296 >      enc[j].resize(nMax);
1297 >      els[j].resize(nMax);
1298 >      ems[j].resize(nMax);
1299 >      ens[j].resize(nMax);
1300      }
1301      
1302      Vector3d t( 2.0 * M_PI );
1303      t.Vdiv(t, box);
1304  
1258 –    
1305      SimInfo::MoleculeIterator mi;
1306      Molecule::AtomIterator ai;
1307      int i;
1308      Vector3d r;
1309      Vector3d tt;
1264 –    Vector3d w;
1265 –    Vector3d u;
1266 –    Vector3d a;
1267 –    Vector3d b;
1310      
1311      for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1312           mol = info_->nextMolecule(mi)) {
# Line 1277 | Line 1319 | namespace OpenMD {
1319          
1320          tt.Vmul(t, r);
1321  
1322 <        
1323 <        eCos[1][i] = Vector3d(1.0, 1.0, 1.0);
1324 <        eSin[1][i] = Vector3d(0.0, 0.0, 0.0);
1325 <        eCos[2][i] = Vector3d(cos(tt.x()), cos(tt.y()), cos(tt.z()));
1326 <        eSin[2][i] = Vector3d(sin(tt.x()), sin(tt.y()), sin(tt.z()));
1322 >        elc[1][i] = 1.0;
1323 >        emc[1][i] = 1.0;
1324 >        enc[1][i] = 1.0;
1325 >        els[1][i] = 0.0;
1326 >        ems[1][i] = 0.0;
1327 >        ens[1][i] = 0.0;
1328  
1329 <        u = eCos[2][i];
1330 <        w = eSin[2][i];
1329 >        elc[2][i] = cos(tt.x());
1330 >        emc[2][i] = cos(tt.y());
1331 >        enc[2][i] = cos(tt.z());
1332 >        els[2][i] = sin(tt.x());
1333 >        ems[2][i] = sin(tt.y());
1334 >        ens[2][i] = sin(tt.z());
1335          
1336          for(int l = 3; l <= kLimit; l++) {
1337 <          eCos[l][i].x() = eCos[l-1][i].x()*eCos[2][i].x() - eSin[l-1][i].x()*eSin[2][i].x();
1338 <          eCos[l][i].y() = eCos[l-1][i].y()*eCos[2][i].y() - eSin[l-1][i].y()*eSin[2][i].y();
1339 <          eCos[l][i].z() = eCos[l-1][i].z()*eCos[2][i].z() - eSin[l-1][i].z()*eSin[2][i].z();
1340 <          
1341 <          eSin[l][i].x() = eSin[l-1][i].x()*eCos[2][i].x() + eCos[l-1][i].x()*eSin[2][i].x();
1342 <          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 <        
1337 >          elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i];
1338 >          emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i];
1339 >          enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i];
1340 >          els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i];
1341 >          ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i];
1342 >          ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i];
1343          }
1344        }
1345      }
# Line 1346 | Line 1383 | namespace OpenMD {
1383      std::vector<RealType> qks(nMax, 0.0);
1384      std::vector<Vector3d> dxk(nMax, V3Zero);
1385      std::vector<Vector3d> qxk(nMax, V3Zero);
1386 <    
1386 >    RealType rl, rm, rn;
1387 >    Vector3d kVec;
1388 >    Vector3d Qk;
1389 >    Mat3x3d k2;
1390 >    RealType ckcs, ckss, dkcs, dkss, qkcs, qkss;
1391 >    int atid;
1392 >    ElectrostaticAtomData data;
1393 >    RealType C, dk, qk;
1394 >    Vector3d D;
1395 >    Mat3x3d  Q;
1396 >
1397      int mMin = kLimit;
1398      int nMin = kLimit + 1;
1399      for (int l = 1; l <= kLimit; l++) {
1400        int ll = l - 1;
1401 <      RealType rl = xcl * float(ll);
1401 >      rl = xcl * float(ll);
1402        for (int mmm = mMin; mmm <= kLim2; mmm++) {
1403          int mm = mmm - kLimit;
1404          int m = abs(mm) + 1;
1405 <        RealType rm = ycl * float(mm);
1405 >        rm = ycl * float(mm);
1406          // Set temporary products of exponential terms
1407          for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1408               mol = info_->nextMolecule(mi)) {
# Line 1364 | Line 1411 | namespace OpenMD {
1411              
1412              i = atom->getLocalIndex();
1413              if(mm < 0) {
1414 <              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1415 <                     + 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();
1414 >              clm[i]=elc[l][i]*emc[m][i]+els[l][i]*ems[m][i];
1415 >              slm[i]=els[l][i]*emc[m][i]-ems[m][i]*elc[l][i];
1416              } else {
1417 <              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1418 <                     - 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();              
1417 >              clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i];
1418 >              slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i];
1419              }
1420            }
1421          }
1422          for (int nnn = nMin; nnn <= kLim2; nnn++) {
1423            int nn = nnn - kLimit;          
1424            int n = abs(nn) + 1;
1425 <          RealType rn = zcl * float(nn);
1425 >          rn = zcl * float(nn);
1426            // Test on magnitude of k vector:
1427            int kk=ll*ll + mm*mm + nn*nn;
1428            if(kk <= kSqLim) {
1429 <            Vector3d kVec = Vector3d(rl, rm, rn);
1430 <            Mat3x3d  k2 = outProduct(kVec, kVec);
1429 >            kVec = Vector3d(rl, rm, rn);
1430 >            k2 = outProduct(kVec, kVec);
1431              // Calculate exp(ikr) terms
1432              for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1433                   mol = info_->nextMolecule(mi)) {
# Line 1393 | Line 1436 | namespace OpenMD {
1436                  i = atom->getLocalIndex();
1437                  
1438                  if (nn < 0) {
1439 <                  ckr[i]=clm[i]*eCos[n][i].z()+slm[i]*eSin[n][i].z();
1440 <                  skr[i]=slm[i]*eCos[n][i].z()-clm[i]*eSin[n][i].z();
1439 >                  ckr[i]=clm[i]*enc[n][i]+slm[i]*ens[n][i];
1440 >                  skr[i]=slm[i]*enc[n][i]-clm[i]*ens[n][i];
1441 >
1442                  } else {
1443 <                  ckr[i]=clm[i]*eCos[n][i].z()-slm[i]*eSin[n][i].z();
1444 <                  skr[i]=slm[i]*eCos[n][i].z()+clm[i]*eSin[n][i].z();
1443 >                  ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i];
1444 >                  skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i];
1445                  }
1446                }
1447              }
# Line 1410 | Line 1454 | namespace OpenMD {
1454                    atom = mol->nextAtom(ai)) {
1455                  i = atom->getLocalIndex();
1456                  int atid = atom->getAtomType()->getIdent();
1457 <                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1457 >                data = ElectrostaticMap[Etids[atid]];
1458                                
1459                  if (data.is_Charge) {
1460 <                  RealType C = data.fixedCharge;
1460 >                  C = data.fixedCharge;
1461                    if (atom->isFluctuatingCharge()) C += atom->getFlucQPos();
1462                    ckc[i] = C * ckr[i];
1463                    cks[i] = C * skr[i];
1464                  }
1465                  
1466                  if (data.is_Dipole) {
1467 <                  Vector3d D = atom->getDipole() * mPoleConverter;
1468 <                  RealType dk = dot(D, kVec);
1467 >                  D = atom->getDipole() * mPoleConverter;
1468 >                  dk = dot(D, kVec);
1469                    dxk[i] = cross(D, kVec);
1470                    dkc[i] = dk * ckr[i];
1471                    dks[i] = dk * skr[i];
1472                  }
1473                  if (data.is_Quadrupole) {
1474 <                  Mat3x3d Q = atom->getQuadrupole();
1475 <                  Q *= mPoleConverter;
1476 <                  RealType qk = - doubleDot(Q, k2);
1477 <                  // RealType qk = -( Q * k2 ).trace();
1434 <                  qxk[i] = -2.0 * cross(k2, Q);
1474 >                  Q = atom->getQuadrupole() * mPoleConverter;
1475 >                  Qk = Q * kVec;                  
1476 >                  qk = dot(kVec, Qk);
1477 >                  qxk[i] = -cross(kVec, Qk);
1478                    qkc[i] = qk * ckr[i];
1479                    qks[i] = qk * skr[i];
1480                  }              
# Line 1440 | Line 1483 | namespace OpenMD {
1483  
1484              // calculate vector sums
1485              
1486 <            RealType ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1487 <            RealType ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1488 <            RealType dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1489 <            RealType dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1490 <            RealType qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1491 <            RealType qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1449 <
1486 >            ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1487 >            ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1488 >            dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1489 >            dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1490 >            qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1491 >            qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1492              
1493   #ifdef IS_MPI
1494 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE,
1495 <                                      MPI::SUM);
1496 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckss, 1, MPI::REALTYPE,
1497 <                                      MPI::SUM);
1498 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkcs, 1, MPI::REALTYPE,
1499 <                                      MPI::SUM);
1500 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkss, 1, MPI::REALTYPE,
1501 <                                      MPI::SUM);
1502 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkcs, 1, MPI::REALTYPE,
1503 <                                      MPI::SUM);
1504 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkss, 1, MPI::REALTYPE,
1505 <                                      MPI::SUM);
1494 >            MPI_Allreduce(MPI_IN_PLACE, &ckcs, 1, MPI_REALTYPE,
1495 >                          MPI_SUM, MPI_COMM_WORLD);
1496 >            MPI_Allreduce(MPI_IN_PLACE, &ckss, 1, MPI_REALTYPE,
1497 >                          MPI_SUM, MPI_COMM_WORLD);
1498 >            MPI_Allreduce(MPI_IN_PLACE, &dkcs, 1, MPI_REALTYPE,
1499 >                          MPI_SUM, MPI_COMM_WORLD);
1500 >            MPI_Allreduce(MPI_IN_PLACE, &dkss, 1, MPI_REALTYPE,
1501 >                          MPI_SUM, MPI_COMM_WORLD);
1502 >            MPI_Allreduce(MPI_IN_PLACE, &qkcs, 1, MPI_REALTYPE,
1503 >                          MPI_SUM, MPI_COMM_WORLD);
1504 >            MPI_Allreduce(MPI_IN_PLACE, &qkss, 1, MPI_REALTYPE,
1505 >                          MPI_SUM, MPI_COMM_WORLD);
1506   #endif        
1507              
1508              // Accumulate potential energy and virial contribution:
1509  
1510 <            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1511 <                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkss));
1510 >            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1511 >                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs));
1512  
1513 <            kVir -= 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss
1514 <                                         +4.0*(ckss*dkcs-ckcs*dkss)
1515 <                                         +3.0*(dkcs*dkcs+dkss*dkss)
1516 <                                         -6.0*(ckss*qkss+ckcs*qkcs)
1517 <                                         +8.0*(dkss*qkcs-dkcs*qkss)
1518 <                                         +5.0*(qkss*qkss+qkcs*qkcs));
1513 >            kVir += 2.0 * rvol  * AK[kk]*(ckcs*ckcs+ckss*ckss
1514 >                                          +4.0*(ckss*dkcs-ckcs*dkss)
1515 >                                          +3.0*(dkcs*dkcs+dkss*dkss)
1516 >                                          -6.0*(ckss*qkss+ckcs*qkcs)
1517 >                                          +8.0*(dkss*qkcs-dkcs*qkss)
1518 >                                          +5.0*(qkss*qkss+qkcs*qkcs));
1519              
1520              // Calculate force and torque for each site:
1521              
# Line 1483 | Line 1525 | namespace OpenMD {
1525                    atom = mol->nextAtom(ai)) {
1526                  
1527                  i = atom->getLocalIndex();
1528 <                int atid = atom->getAtomType()->getIdent();
1529 <                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1530 <                
1528 >                atid = atom->getAtomType()->getIdent();
1529 >                data = ElectrostaticMap[Etids[atid]];
1530 >
1531                  RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs)
1532                                       - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss));
1533                  RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
1534                                           -ckr[i]*(ckss+dkcs-qkss));
1535 <                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)+
1536 <                                             skr[i]*(ckss+dkcs-qkss));
1535 >                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
1536 >                                            +skr[i]*(ckss+dkcs-qkss));
1537                
1538                  atom->addFrc( 4.0 * rvol * qfrc * kVec );
1539 <                
1539 >
1540 >                if (atom->isFluctuatingCharge()) {
1541 >                  atom->addFlucQFrc( - 2.0 * rvol * qtrq2 );
1542 >                }
1543 >                  
1544                  if (data.is_Dipole) {
1545                    atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] );
1546                  }
# Line 1509 | Line 1555 | namespace OpenMD {
1555        }
1556        mMin = 1;
1557      }
1558 <    cerr << "kPot = " << kPot << "\n";
1513 <    pot[ELECTROSTATIC_FAMILY] += kPot;  
1558 >    pot += kPot;  
1559    }
1560 +
1561 +  void Electrostatic::getSitePotentials(Atom* a1, Atom* a2, bool excluded,
1562 +                                        RealType &spot1, RealType &spot2) {
1563 +
1564 +    if (!initialized_) {
1565 +      cerr << "initializing\n";
1566 +      initialize();
1567 +      cerr << "done\n";
1568 +    }
1569 +
1570 +    const RealType mPoleConverter = 0.20819434;
1571 +
1572 +    AtomType* atype1 = a1->getAtomType();
1573 +    AtomType* atype2 = a2->getAtomType();
1574 +    int atid1 = atype1->getIdent();
1575 +    int atid2 = atype2->getIdent();
1576 +    data1 = ElectrostaticMap[Etids[atid1]];
1577 +    data2 = ElectrostaticMap[Etids[atid2]];
1578 +
1579 +    Pa = 0.0;  // Site potential at site a
1580 +    Pb = 0.0;  // Site potential at site b
1581 +
1582 +    Vector3d d = a2->getPos() - a1->getPos();
1583 +    info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(d);
1584 +    RealType rij = d.length();
1585 +    // some variables we'll need independent of electrostatic type:
1586 +
1587 +    RealType ri = 1.0 /  rij;
1588 +    rhat = d  * ri;
1589 +      
1590 +
1591 +    if ((rij >= cutoffRadius_) || excluded) {
1592 +      spot1 = 0.0;
1593 +      spot2 = 0.0;
1594 +      return;
1595 +    }
1596 +
1597 +    // logicals
1598 +
1599 +    a_is_Charge = data1.is_Charge;
1600 +    a_is_Dipole = data1.is_Dipole;
1601 +    a_is_Quadrupole = data1.is_Quadrupole;
1602 +    a_is_Fluctuating = data1.is_Fluctuating;
1603 +
1604 +    b_is_Charge = data2.is_Charge;
1605 +    b_is_Dipole = data2.is_Dipole;
1606 +    b_is_Quadrupole = data2.is_Quadrupole;
1607 +    b_is_Fluctuating = data2.is_Fluctuating;
1608 +
1609 +    // Obtain all of the required radial function values from the
1610 +    // spline structures:
1611 +    
1612 +
1613 +    if (a_is_Charge || b_is_Charge) {
1614 +      v01 = v01s->getValueAt(rij);
1615 +    }
1616 +    if (a_is_Dipole || b_is_Dipole) {
1617 +      v11 = v11s->getValueAt(rij);
1618 +      v11or = ri * v11;
1619 +    }
1620 +    if (a_is_Quadrupole || b_is_Quadrupole) {
1621 +      v21 = v21s->getValueAt(rij);
1622 +      v22 = v22s->getValueAt(rij);
1623 +      v22or = ri * v22;
1624 +    }      
1625 +
1626 +    if (a_is_Charge) {
1627 +      C_a = data1.fixedCharge;
1628 +      
1629 +      if (a_is_Fluctuating) {
1630 +        C_a += a1->getFlucQPos();
1631 +      }
1632 +      
1633 +      Pb += C_a *  pre11_ * v01;      
1634 +    }
1635 +    
1636 +    if (a_is_Dipole) {
1637 +      D_a = a1->getDipole() * mPoleConverter;
1638 +      rdDa = dot(rhat, D_a);
1639 +      Pb +=  pre12_ * v11 * rdDa;      
1640 +    }
1641 +    
1642 +    if (a_is_Quadrupole) {
1643 +      Q_a = a1->getQuadrupole() * mPoleConverter;
1644 +      trQa =  Q_a.trace();
1645 +      Qar =   Q_a * rhat;
1646 +      rdQar = dot(rhat, Qar);
1647 +      Pb += pre14_ * (v21 * trQa + v22 * rdQar);      
1648 +    }
1649 +    
1650 +    if (b_is_Charge) {
1651 +      C_b = data2.fixedCharge;
1652 +      
1653 +      if (b_is_Fluctuating)
1654 +        C_b += a2->getFlucQPos();
1655 +      
1656 +      Pa += C_b *  pre11_ * v01;
1657 +    }
1658 +    
1659 +    if (b_is_Dipole) {
1660 +      D_b = a2->getDipole() * mPoleConverter;
1661 +      rdDb = dot(rhat, D_b);
1662 +      Pa += pre12_ * v11 * rdDb;
1663 +    }
1664 +    
1665 +    if (b_is_Quadrupole) {
1666 +      Q_a = a2->getQuadrupole() * mPoleConverter;
1667 +      trQb =  Q_b.trace();
1668 +      Qbr =   Q_b * rhat;
1669 +      rdQbr = dot(rhat, Qbr);
1670 +      Pa += pre14_ * (v21 * trQb + v22 * rdQbr);      
1671 +    }
1672 +
1673 +    spot1 = Pa;
1674 +    spot2 = Pb;
1675 +  }
1676   }

Diff Legend

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