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 1929 by gezelter, Mon Aug 19 13:12:00 2013 UTC

# Line 57 | Line 57
57   #include "math/erfc.hpp"
58   #include "math/SquareMatrix.hpp"
59   #include "primitives/Molecule.hpp"
60 <
60 > #ifdef IS_MPI
61 > #include <mpi.h>
62 > #endif
63  
64   namespace OpenMD {
65    
# Line 886 | Line 888 | namespace OpenMD {
888                          + rdQbr * rhat * (dv22 - 2.0*v22or));
889      }
890      
891 <    if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
891 >    if ((a_is_Fluctuating || b_is_Fluctuating)
892 >        && (idat.excluded || idat.sameRegion)) {
893        J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
894      }    
895      
# Line 907 | Line 910 | namespace OpenMD {
910          }
911          
912          // Fluctuating charge forces are handled via Coulomb integrals
913 <        // for excluded pairs (i.e. those connected via bonds) and
914 <        // with the standard charge-charge interaction otherwise.
913 >        // for excluded pairs (i.e. those connected via bonds), atoms
914 >        // within the same region (for metals) and with the standard
915 >        // charge-charge interaction otherwise.
916  
917 <        if (idat.excluded) {          
917 >        if (idat.excluded || idat.sameRegion) {          
918            if (a_is_Fluctuating || b_is_Fluctuating) {
919              coulInt = J->getValueAt( *(idat.rij) );
920              if (a_is_Fluctuating)  dUdCa += coulInt * C_b;
921              if (b_is_Fluctuating)  dUdCb += coulInt * C_a;
922 <            excluded_Pot += C_a * C_b * coulInt;
922 >            if (idat.excluded)
923 >              excluded_Pot += C_a * C_b * coulInt;
924            }          
925          } else {
926            if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
927 <          if (a_is_Fluctuating) dUdCb += C_a * pref * v01;
927 >          if (b_is_Fluctuating) dUdCb += C_a * pref * v01;
928          }
929        }
930  
# Line 1140 | Line 1145 | namespace OpenMD {
1145          
1146      if (i_is_Fluctuating) {
1147        C_a += *(sdat.flucQ);
1148 <      // dVdFQ is really a force, so this is negative the derivative
1144 <      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1148 >      *(sdat.flucQfrc) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1149        (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1150          (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1151      }
# Line 1194 | Line 1198 | namespace OpenMD {
1198    }
1199  
1200  
1201 <  void Electrostatic::ReciprocalSpaceSum(potVec& pot) {
1201 >  void Electrostatic::ReciprocalSpaceSum(RealType& pot) {
1202      
1203      RealType kPot = 0.0;
1204      RealType kVir = 0.0;
# Line 1240 | Line 1244 | namespace OpenMD {
1244      
1245      // Calculate and store exponential factors  
1246      
1247 <    vector<vector<Vector3d> > eCos;
1248 <    vector<vector<Vector3d> > eSin;
1247 >    vector<vector<RealType> > elc;
1248 >    vector<vector<RealType> > emc;
1249 >    vector<vector<RealType> > enc;
1250 >    vector<vector<RealType> > els;
1251 >    vector<vector<RealType> > ems;
1252 >    vector<vector<RealType> > ens;
1253 >
1254      
1255      int nMax = info_->getNAtoms();
1256      
1257 <    eCos.resize(kLimit+1);
1258 <    eSin.resize(kLimit+1);
1257 >    elc.resize(kLimit+1);
1258 >    emc.resize(kLimit+1);
1259 >    enc.resize(kLimit+1);
1260 >    els.resize(kLimit+1);
1261 >    ems.resize(kLimit+1);
1262 >    ens.resize(kLimit+1);
1263 >
1264      for (int j = 0; j < kLimit+1; j++) {
1265 <      eCos[j].resize(nMax);
1266 <      eSin[j].resize(nMax);
1265 >      elc[j].resize(nMax);
1266 >      emc[j].resize(nMax);
1267 >      enc[j].resize(nMax);
1268 >      els[j].resize(nMax);
1269 >      ems[j].resize(nMax);
1270 >      ens[j].resize(nMax);
1271      }
1272      
1273      Vector3d t( 2.0 * M_PI );
# Line 1261 | Line 1279 | namespace OpenMD {
1279      int i;
1280      Vector3d r;
1281      Vector3d tt;
1264 –    Vector3d w;
1265 –    Vector3d u;
1266 –    Vector3d a;
1267 –    Vector3d b;
1282      
1283      for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1284           mol = info_->nextMolecule(mi)) {
# Line 1277 | Line 1291 | namespace OpenMD {
1291          
1292          tt.Vmul(t, r);
1293  
1294 <        
1295 <        eCos[1][i] = Vector3d(1.0, 1.0, 1.0);
1296 <        eSin[1][i] = Vector3d(0.0, 0.0, 0.0);
1297 <        eCos[2][i] = Vector3d(cos(tt.x()), cos(tt.y()), cos(tt.z()));
1298 <        eSin[2][i] = Vector3d(sin(tt.x()), sin(tt.y()), sin(tt.z()));
1294 >        elc[1][i] = 1.0;
1295 >        emc[1][i] = 1.0;
1296 >        enc[1][i] = 1.0;
1297 >        els[1][i] = 0.0;
1298 >        ems[1][i] = 0.0;
1299 >        ens[1][i] = 0.0;
1300  
1301 <        u = eCos[2][i];
1302 <        w = eSin[2][i];
1301 >        elc[2][i] = cos(tt.x());
1302 >        emc[2][i] = cos(tt.y());
1303 >        enc[2][i] = cos(tt.z());
1304 >        els[2][i] = sin(tt.x());
1305 >        ems[2][i] = sin(tt.y());
1306 >        ens[2][i] = sin(tt.z());
1307          
1308          for(int l = 3; l <= kLimit; l++) {
1309 <          eCos[l][i].x() = eCos[l-1][i].x()*eCos[2][i].x() - eSin[l-1][i].x()*eSin[2][i].x();
1310 <          eCos[l][i].y() = eCos[l-1][i].y()*eCos[2][i].y() - eSin[l-1][i].y()*eSin[2][i].y();
1311 <          eCos[l][i].z() = eCos[l-1][i].z()*eCos[2][i].z() - eSin[l-1][i].z()*eSin[2][i].z();
1312 <          
1313 <          eSin[l][i].x() = eSin[l-1][i].x()*eCos[2][i].x() + eCos[l-1][i].x()*eSin[2][i].x();
1314 <          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 <        
1309 >          elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i];
1310 >          emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i];
1311 >          enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i];
1312 >          els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i];
1313 >          ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i];
1314 >          ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i];
1315          }
1316        }
1317      }
# Line 1346 | Line 1355 | namespace OpenMD {
1355      std::vector<RealType> qks(nMax, 0.0);
1356      std::vector<Vector3d> dxk(nMax, V3Zero);
1357      std::vector<Vector3d> qxk(nMax, V3Zero);
1358 <    
1358 >    RealType rl, rm, rn;
1359 >    Vector3d kVec;
1360 >    Vector3d Qk;
1361 >    Mat3x3d k2;
1362 >    RealType ckcs, ckss, dkcs, dkss, qkcs, qkss;
1363 >    int atid;
1364 >    ElectrostaticAtomData data;
1365 >    RealType C, dk, qk;
1366 >    Vector3d D;
1367 >    Mat3x3d  Q;
1368 >
1369      int mMin = kLimit;
1370      int nMin = kLimit + 1;
1371      for (int l = 1; l <= kLimit; l++) {
1372        int ll = l - 1;
1373 <      RealType rl = xcl * float(ll);
1373 >      rl = xcl * float(ll);
1374        for (int mmm = mMin; mmm <= kLim2; mmm++) {
1375          int mm = mmm - kLimit;
1376          int m = abs(mm) + 1;
1377 <        RealType rm = ycl * float(mm);
1377 >        rm = ycl * float(mm);
1378          // Set temporary products of exponential terms
1379          for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1380               mol = info_->nextMolecule(mi)) {
# Line 1364 | Line 1383 | namespace OpenMD {
1383              
1384              i = atom->getLocalIndex();
1385              if(mm < 0) {
1386 <              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1387 <                     + 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();
1386 >              clm[i]=elc[l][i]*emc[m][i]+els[l][i]*ems[m][i];
1387 >              slm[i]=els[l][i]*emc[m][i]-ems[m][i]*elc[l][i];
1388              } else {
1389 <              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1390 <                     - 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();              
1389 >              clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i];
1390 >              slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i];
1391              }
1392            }
1393          }
1394          for (int nnn = nMin; nnn <= kLim2; nnn++) {
1395            int nn = nnn - kLimit;          
1396            int n = abs(nn) + 1;
1397 <          RealType rn = zcl * float(nn);
1397 >          rn = zcl * float(nn);
1398            // Test on magnitude of k vector:
1399            int kk=ll*ll + mm*mm + nn*nn;
1400            if(kk <= kSqLim) {
1401 <            Vector3d kVec = Vector3d(rl, rm, rn);
1402 <            Mat3x3d  k2 = outProduct(kVec, kVec);
1401 >            kVec = Vector3d(rl, rm, rn);
1402 >            k2 = outProduct(kVec, kVec);
1403              // Calculate exp(ikr) terms
1404              for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1405                   mol = info_->nextMolecule(mi)) {
# Line 1393 | Line 1408 | namespace OpenMD {
1408                  i = atom->getLocalIndex();
1409                  
1410                  if (nn < 0) {
1411 <                  ckr[i]=clm[i]*eCos[n][i].z()+slm[i]*eSin[n][i].z();
1412 <                  skr[i]=slm[i]*eCos[n][i].z()-clm[i]*eSin[n][i].z();
1411 >                  ckr[i]=clm[i]*enc[n][i]+slm[i]*ens[n][i];
1412 >                  skr[i]=slm[i]*enc[n][i]-clm[i]*ens[n][i];
1413 >
1414                  } else {
1415 <                  ckr[i]=clm[i]*eCos[n][i].z()-slm[i]*eSin[n][i].z();
1416 <                  skr[i]=slm[i]*eCos[n][i].z()+clm[i]*eSin[n][i].z();
1415 >                  ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i];
1416 >                  skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i];
1417                  }
1418                }
1419              }
# Line 1410 | Line 1426 | namespace OpenMD {
1426                    atom = mol->nextAtom(ai)) {
1427                  i = atom->getLocalIndex();
1428                  int atid = atom->getAtomType()->getIdent();
1429 <                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1429 >                data = ElectrostaticMap[Etids[atid]];
1430                                
1431                  if (data.is_Charge) {
1432 <                  RealType C = data.fixedCharge;
1432 >                  C = data.fixedCharge;
1433                    if (atom->isFluctuatingCharge()) C += atom->getFlucQPos();
1434                    ckc[i] = C * ckr[i];
1435                    cks[i] = C * skr[i];
1436                  }
1437                  
1438                  if (data.is_Dipole) {
1439 <                  Vector3d D = atom->getDipole() * mPoleConverter;
1440 <                  RealType dk = dot(D, kVec);
1439 >                  D = atom->getDipole() * mPoleConverter;
1440 >                  dk = dot(D, kVec);
1441                    dxk[i] = cross(D, kVec);
1442                    dkc[i] = dk * ckr[i];
1443                    dks[i] = dk * skr[i];
1444                  }
1445                  if (data.is_Quadrupole) {
1446 <                  Mat3x3d Q = atom->getQuadrupole();
1446 >                  Q = atom->getQuadrupole();
1447                    Q *= mPoleConverter;
1448 <                  RealType qk = - doubleDot(Q, k2);
1449 <                  // RealType qk = -( Q * k2 ).trace();
1450 <                  qxk[i] = -2.0 * cross(k2, Q);
1448 >                  Qk = Q * kVec;
1449 >                  qk = dot(kVec, Qk);
1450 >                  qxk[i] = cross(kVec, Qk);
1451                    qkc[i] = qk * ckr[i];
1452                    qks[i] = qk * skr[i];
1453                  }              
# Line 1440 | Line 1456 | namespace OpenMD {
1456  
1457              // calculate vector sums
1458              
1459 <            RealType ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1460 <            RealType ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1461 <            RealType dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1462 <            RealType dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1463 <            RealType qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1464 <            RealType qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1449 <
1459 >            ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1460 >            ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1461 >            dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1462 >            dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1463 >            qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1464 >            qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1465              
1466   #ifdef IS_MPI
1467              MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE,
# Line 1465 | Line 1480 | namespace OpenMD {
1480              
1481              // Accumulate potential energy and virial contribution:
1482  
1483 <            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1484 <                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkss));
1483 >            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1484 >                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs));
1485  
1486 <            kVir -= 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss
1487 <                                         +4.0*(ckss*dkcs-ckcs*dkss)
1488 <                                         +3.0*(dkcs*dkcs+dkss*dkss)
1489 <                                         -6.0*(ckss*qkss+ckcs*qkcs)
1490 <                                         +8.0*(dkss*qkcs-dkcs*qkss)
1491 <                                         +5.0*(qkss*qkss+qkcs*qkcs));
1486 >            kVir += 2.0 * rvol  * AK[kk]*(ckcs*ckcs+ckss*ckss
1487 >                                          +4.0*(ckss*dkcs-ckcs*dkss)
1488 >                                          +3.0*(dkcs*dkcs+dkss*dkss)
1489 >                                          -6.0*(ckss*qkss+ckcs*qkcs)
1490 >                                          +8.0*(dkss*qkcs-dkcs*qkss)
1491 >                                          +5.0*(qkss*qkss+qkcs*qkcs));
1492              
1493              // Calculate force and torque for each site:
1494              
# Line 1483 | Line 1498 | namespace OpenMD {
1498                    atom = mol->nextAtom(ai)) {
1499                  
1500                  i = atom->getLocalIndex();
1501 <                int atid = atom->getAtomType()->getIdent();
1502 <                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1503 <                
1501 >                atid = atom->getAtomType()->getIdent();
1502 >                data = ElectrostaticMap[Etids[atid]];
1503 >
1504                  RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs)
1505                                       - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss));
1506                  RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
1507                                           -ckr[i]*(ckss+dkcs-qkss));
1508 <                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)+
1509 <                                             skr[i]*(ckss+dkcs-qkss));
1508 >                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
1509 >                                             +skr[i]*(ckss+dkcs-qkss));
1510                
1511                  atom->addFrc( 4.0 * rvol * qfrc * kVec );
1512                  
# Line 1509 | Line 1524 | namespace OpenMD {
1524        }
1525        mMin = 1;
1526      }
1527 <    cerr << "kPot = " << kPot << "\n";
1513 <    pot[ELECTROSTATIC_FAMILY] += kPot;  
1527 >    pot += kPot;  
1528    }
1529   }

Diff Legend

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