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 1923 by gezelter, Mon Aug 5 16:13:46 2013 UTC vs.
Revision 1933 by gezelter, Fri Aug 23 15:59:23 2013 UTC

# Line 57 | Line 57
57   #include "math/erfc.hpp"
58   #include "math/SquareMatrix.hpp"
59   #include "primitives/Molecule.hpp"
60 + #ifdef IS_MPI
61 + #include <mpi.h>
62 + #endif
63  
61 –
64   namespace OpenMD {
65    
66    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
# Line 764 | Line 766 | namespace OpenMD {
766      // Excluded potential that is still computed for fluctuating charges
767      excluded_Pot= 0.0;
768  
767 –
769      // some variables we'll need independent of electrostatic type:
770  
771      ri = 1.0 /  *(idat.rij);
# Line 885 | 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 905 | 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;
919 <          }          
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 986 | 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);
989 –
990          // Even if we excluded this pair from direct interactions, we
991          // still have the reaction-field-mediated dipole-dipole
992          // interaction:
# Line 1046 | 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 1077 | Line 1077 | namespace OpenMD {
1077          //             + 4.0 * cross(rhat, QbQar)
1078  
1079          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1080 –
1080        }
1081      }
1082  
# Line 1140 | 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
1144 <      *(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 1194 | Line 1192 | namespace OpenMD {
1192    }
1193  
1194  
1195 <  void Electrostatic::ReciprocalSpaceSum(potVec& pot) {
1195 >  void Electrostatic::ReciprocalSpaceSum(RealType& pot) {
1196      
1197      RealType kPot = 0.0;
1198      RealType kVir = 0.0;
# Line 1240 | Line 1238 | namespace OpenMD {
1238      
1239      // Calculate and store exponential factors  
1240      
1241 <    vector<vector<Vector3d> > eCos;
1242 <    vector<vector<Vector3d> > eSin;
1241 >    vector<vector<RealType> > elc;
1242 >    vector<vector<RealType> > emc;
1243 >    vector<vector<RealType> > enc;
1244 >    vector<vector<RealType> > els;
1245 >    vector<vector<RealType> > ems;
1246 >    vector<vector<RealType> > ens;
1247 >
1248      
1249      int nMax = info_->getNAtoms();
1250      
1251 <    eCos.resize(kLimit+1);
1252 <    eSin.resize(kLimit+1);
1251 >    elc.resize(kLimit+1);
1252 >    emc.resize(kLimit+1);
1253 >    enc.resize(kLimit+1);
1254 >    els.resize(kLimit+1);
1255 >    ems.resize(kLimit+1);
1256 >    ens.resize(kLimit+1);
1257 >
1258      for (int j = 0; j < kLimit+1; j++) {
1259 <      eCos[j].resize(nMax);
1260 <      eSin[j].resize(nMax);
1259 >      elc[j].resize(nMax);
1260 >      emc[j].resize(nMax);
1261 >      enc[j].resize(nMax);
1262 >      els[j].resize(nMax);
1263 >      ems[j].resize(nMax);
1264 >      ens[j].resize(nMax);
1265      }
1266      
1267      Vector3d t( 2.0 * M_PI );
# Line 1261 | Line 1273 | namespace OpenMD {
1273      int i;
1274      Vector3d r;
1275      Vector3d tt;
1264 –    Vector3d w;
1265 –    Vector3d u;
1266 –    Vector3d a;
1267 –    Vector3d b;
1276      
1277      for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1278           mol = info_->nextMolecule(mi)) {
# Line 1277 | Line 1285 | namespace OpenMD {
1285          
1286          tt.Vmul(t, r);
1287  
1288 <        
1289 <        eCos[1][i] = Vector3d(1.0, 1.0, 1.0);
1290 <        eSin[1][i] = Vector3d(0.0, 0.0, 0.0);
1291 <        eCos[2][i] = Vector3d(cos(tt.x()), cos(tt.y()), cos(tt.z()));
1292 <        eSin[2][i] = Vector3d(sin(tt.x()), sin(tt.y()), sin(tt.z()));
1288 >        elc[1][i] = 1.0;
1289 >        emc[1][i] = 1.0;
1290 >        enc[1][i] = 1.0;
1291 >        els[1][i] = 0.0;
1292 >        ems[1][i] = 0.0;
1293 >        ens[1][i] = 0.0;
1294  
1295 <        u = eCos[2][i];
1296 <        w = eSin[2][i];
1295 >        elc[2][i] = cos(tt.x());
1296 >        emc[2][i] = cos(tt.y());
1297 >        enc[2][i] = cos(tt.z());
1298 >        els[2][i] = sin(tt.x());
1299 >        ems[2][i] = sin(tt.y());
1300 >        ens[2][i] = sin(tt.z());
1301          
1302          for(int l = 3; l <= kLimit; l++) {
1303 <          eCos[l][i].x() = eCos[l-1][i].x()*eCos[2][i].x() - eSin[l-1][i].x()*eSin[2][i].x();
1304 <          eCos[l][i].y() = eCos[l-1][i].y()*eCos[2][i].y() - eSin[l-1][i].y()*eSin[2][i].y();
1305 <          eCos[l][i].z() = eCos[l-1][i].z()*eCos[2][i].z() - eSin[l-1][i].z()*eSin[2][i].z();
1306 <          
1307 <          eSin[l][i].x() = eSin[l-1][i].x()*eCos[2][i].x() + eCos[l-1][i].x()*eSin[2][i].x();
1308 <          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 <        
1303 >          elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i];
1304 >          emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i];
1305 >          enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i];
1306 >          els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i];
1307 >          ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i];
1308 >          ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i];
1309          }
1310        }
1311      }
# Line 1346 | Line 1349 | namespace OpenMD {
1349      std::vector<RealType> qks(nMax, 0.0);
1350      std::vector<Vector3d> dxk(nMax, V3Zero);
1351      std::vector<Vector3d> qxk(nMax, V3Zero);
1352 <    
1352 >    RealType rl, rm, rn;
1353 >    Vector3d kVec;
1354 >    Vector3d Qk;
1355 >    Mat3x3d k2;
1356 >    RealType ckcs, ckss, dkcs, dkss, qkcs, qkss;
1357 >    int atid;
1358 >    ElectrostaticAtomData data;
1359 >    RealType C, dk, qk;
1360 >    Vector3d D;
1361 >    Mat3x3d  Q;
1362 >
1363      int mMin = kLimit;
1364      int nMin = kLimit + 1;
1365      for (int l = 1; l <= kLimit; l++) {
1366        int ll = l - 1;
1367 <      RealType rl = xcl * float(ll);
1367 >      rl = xcl * float(ll);
1368        for (int mmm = mMin; mmm <= kLim2; mmm++) {
1369          int mm = mmm - kLimit;
1370          int m = abs(mm) + 1;
1371 <        RealType rm = ycl * float(mm);
1371 >        rm = ycl * float(mm);
1372          // Set temporary products of exponential terms
1373          for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1374               mol = info_->nextMolecule(mi)) {
# Line 1364 | Line 1377 | namespace OpenMD {
1377              
1378              i = atom->getLocalIndex();
1379              if(mm < 0) {
1380 <              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1381 <                     + 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();
1380 >              clm[i]=elc[l][i]*emc[m][i]+els[l][i]*ems[m][i];
1381 >              slm[i]=els[l][i]*emc[m][i]-ems[m][i]*elc[l][i];
1382              } else {
1383 <              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1384 <                     - 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();              
1383 >              clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i];
1384 >              slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i];
1385              }
1386            }
1387          }
1388          for (int nnn = nMin; nnn <= kLim2; nnn++) {
1389            int nn = nnn - kLimit;          
1390            int n = abs(nn) + 1;
1391 <          RealType rn = zcl * float(nn);
1391 >          rn = zcl * float(nn);
1392            // Test on magnitude of k vector:
1393            int kk=ll*ll + mm*mm + nn*nn;
1394            if(kk <= kSqLim) {
1395 <            Vector3d kVec = Vector3d(rl, rm, rn);
1396 <            Mat3x3d  k2 = outProduct(kVec, kVec);
1395 >            kVec = Vector3d(rl, rm, rn);
1396 >            k2 = outProduct(kVec, kVec);
1397              // Calculate exp(ikr) terms
1398              for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1399                   mol = info_->nextMolecule(mi)) {
# Line 1393 | Line 1402 | namespace OpenMD {
1402                  i = atom->getLocalIndex();
1403                  
1404                  if (nn < 0) {
1405 <                  ckr[i]=clm[i]*eCos[n][i].z()+slm[i]*eSin[n][i].z();
1406 <                  skr[i]=slm[i]*eCos[n][i].z()-clm[i]*eSin[n][i].z();
1405 >                  ckr[i]=clm[i]*enc[n][i]+slm[i]*ens[n][i];
1406 >                  skr[i]=slm[i]*enc[n][i]-clm[i]*ens[n][i];
1407 >
1408                  } else {
1409 <                  ckr[i]=clm[i]*eCos[n][i].z()-slm[i]*eSin[n][i].z();
1410 <                  skr[i]=slm[i]*eCos[n][i].z()+clm[i]*eSin[n][i].z();
1409 >                  ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i];
1410 >                  skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i];
1411                  }
1412                }
1413              }
# Line 1410 | Line 1420 | namespace OpenMD {
1420                    atom = mol->nextAtom(ai)) {
1421                  i = atom->getLocalIndex();
1422                  int atid = atom->getAtomType()->getIdent();
1423 <                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1423 >                data = ElectrostaticMap[Etids[atid]];
1424                                
1425                  if (data.is_Charge) {
1426 <                  RealType C = data.fixedCharge;
1426 >                  C = data.fixedCharge;
1427                    if (atom->isFluctuatingCharge()) C += atom->getFlucQPos();
1428                    ckc[i] = C * ckr[i];
1429                    cks[i] = C * skr[i];
1430                  }
1431                  
1432                  if (data.is_Dipole) {
1433 <                  Vector3d D = atom->getDipole() * mPoleConverter;
1434 <                  RealType dk = dot(kVec, D);
1435 <                  dxk[i] = cross(kVec, D);
1433 >                  D = atom->getDipole() * mPoleConverter;
1434 >                  dk = dot(D, kVec);
1435 >                  dxk[i] = cross(D, kVec);
1436                    dkc[i] = dk * ckr[i];
1437                    dks[i] = dk * skr[i];
1438                  }
1439                  if (data.is_Quadrupole) {
1440 <                  Mat3x3d Q = atom->getQuadrupole();
1441 <                  Q *= mPoleConverter;
1442 <                  RealType qk = -( Q * k2 ).trace();
1443 <                  qxk[i] = -2.0 * cross(k2, Q);
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];
1445                    qks[i] = qk * skr[i];
1446                  }              
# Line 1439 | Line 1449 | namespace OpenMD {
1449  
1450              // calculate vector sums
1451              
1452 <            RealType ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1453 <            RealType ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1454 <            RealType dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1455 <            RealType dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1456 <            RealType qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1457 <            RealType qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1448 <
1452 >            ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1453 >            ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1454 >            dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1455 >            dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1456 >            qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1457 >            qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1458              
1459   #ifdef IS_MPI
1460              MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE,
# Line 1464 | Line 1473 | namespace OpenMD {
1473              
1474              // Accumulate potential energy and virial contribution:
1475  
1476 <            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1477 <                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkss));
1476 >            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1477 >                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs));
1478  
1479 <            kVir -= 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss
1480 <                                         +4.0*(ckss*dkcs-ckcs*dkss)
1481 <                                         +3.0*(dkcs*dkcs+dkss*dkss)
1482 <                                         -6.0*(ckss*qkss+ckcs*qkcs)
1483 <                                         +8.0*(dkss*qkcs-dkcs*qkss)
1484 <                                         +5.0*(qkss*qkss+qkcs*qkcs));
1479 >            kVir += 2.0 * rvol  * AK[kk]*(ckcs*ckcs+ckss*ckss
1480 >                                          +4.0*(ckss*dkcs-ckcs*dkss)
1481 >                                          +3.0*(dkcs*dkcs+dkss*dkss)
1482 >                                          -6.0*(ckss*qkss+ckcs*qkcs)
1483 >                                          +8.0*(dkss*qkcs-dkcs*qkss)
1484 >                                          +5.0*(qkss*qkss+qkcs*qkcs));
1485              
1486              // Calculate force and torque for each site:
1487              
# Line 1482 | Line 1491 | namespace OpenMD {
1491                    atom = mol->nextAtom(ai)) {
1492                  
1493                  i = atom->getLocalIndex();
1494 <                int atid = atom->getAtomType()->getIdent();
1495 <                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1496 <                
1494 >                atid = atom->getAtomType()->getIdent();
1495 >                data = ElectrostaticMap[Etids[atid]];
1496 >
1497                  RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs)
1498                                       - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss));
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));
1501 >                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
1502 >                                            +skr[i]*(ckss+dkcs-qkss));
1503                
1504                  atom->addFrc( 4.0 * rvol * qfrc * kVec );
1505                  
# Line 1508 | Line 1517 | namespace OpenMD {
1517        }
1518        mMin = 1;
1519      }
1520 <    cerr << "kPot = " << kPot << "\n";
1512 <    pot[ELECTROSTATIC_FAMILY] += kPot;  
1520 >    pot += kPot;  
1521    }
1522   }

Diff Legend

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