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 1915 by gezelter, Mon Jul 29 17:55:17 2013 UTC vs.
Revision 1938 by gezelter, Thu Oct 31 15:32:17 2013 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 58 | Line 62
62   #include "math/SquareMatrix.hpp"
63   #include "primitives/Molecule.hpp"
64  
61 –
65   namespace OpenMD {
66    
67    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
# Line 677 | Line 680 | namespace OpenMD {
680        FQtids[atid] = fqtid;
681        Jij[fqtid].resize(nFlucq_);
682  
683 <      // Now, iterate over all known fluctuating and add to the coulomb integral map:
683 >      // Now, iterate over all known fluctuating and add to the
684 >      // coulomb integral map:
685        
686        std::set<int>::iterator it;
687        for( it = FQtypes.begin(); it != FQtypes.end(); ++it) {    
# Line 763 | Line 767 | namespace OpenMD {
767      // Excluded potential that is still computed for fluctuating charges
768      excluded_Pot= 0.0;
769  
766 –
770      // some variables we'll need independent of electrostatic type:
771  
772      ri = 1.0 /  *(idat.rij);
# Line 884 | Line 887 | namespace OpenMD {
887          Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or
888                          + rdQbr * rhat * (dv22 - 2.0*v22or));
889      }
890 <    
890 >        
891 >
892      if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
893        J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
894      }    
895 <    
895 >
896      if (a_is_Charge) {    
897        
898        if (b_is_Charge) {
899          pref =  pre11_ * *(idat.electroMult);      
900          U  += C_a * C_b * pref * v01;
901          F  += C_a * C_b * pref * dv01 * rhat;
902 <        
902 >
903          // If this is an excluded pair, there are still indirect
904          // interactions via the reaction field we must worry about:
905  
# Line 904 | Line 908 | namespace OpenMD {
908            indirect_Pot += rfContrib;
909            indirect_F   += rfContrib * 2.0 * ri * rhat;
910          }
911 <        
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.
915  
916 <        if (idat.excluded) {          
916 >        if (idat.excluded) {
917            if (a_is_Fluctuating || b_is_Fluctuating) {
918              coulInt = J->getValueAt( *(idat.rij) );
919 <            if (a_is_Fluctuating)  dUdCa += coulInt * C_b;
920 <            if (b_is_Fluctuating)  dUdCb += coulInt * C_a;
921 <            excluded_Pot += C_a * C_b * coulInt;
918 <          }          
919 >            if (a_is_Fluctuating) dUdCa += C_b * coulInt;
920 >            if (b_is_Fluctuating) dUdCb += C_a * coulInt;          
921 >          }
922          } else {
923            if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
924 <          if (a_is_Fluctuating) dUdCb += C_a * pref * v01;
925 <        }
924 >          if (b_is_Fluctuating) dUdCb += C_a * pref * v01;
925 >        }              
926        }
927  
928        if (b_is_Dipole) {
# Line 985 | Line 988 | namespace OpenMD {
988          F  -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat;
989          Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
990          Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
988 –
991          // Even if we excluded this pair from direct interactions, we
992          // still have the reaction-field-mediated dipole-dipole
993          // interaction:
# Line 1045 | Line 1047 | namespace OpenMD {
1047          trQaQb = QaQb.trace();
1048          rQaQb = rhat * QaQb;
1049          QaQbr = QaQb * rhat;
1050 <        QaxQb = cross(Q_a, Q_b);
1050 >        QaxQb = mCross(Q_a, Q_b);
1051          rQaQbr = dot(rQa, Qbr);
1052          rQaxQbr = cross(rQa, Qbr);
1053          
# Line 1076 | Line 1078 | namespace OpenMD {
1078          //             + 4.0 * cross(rhat, QbQar)
1079  
1080          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1079 –
1081        }
1082      }
1083  
# Line 1139 | Line 1140 | namespace OpenMD {
1140          
1141      if (i_is_Fluctuating) {
1142        C_a += *(sdat.flucQ);
1143 <      // dVdFQ is really a force, so this is negative the derivative
1143 <      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1143 >      *(sdat.flucQfrc) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1144        (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1145          (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1146      }
# Line 1165 | Line 1165 | namespace OpenMD {
1165      case esm_SHIFTED_FORCE:
1166      case esm_SHIFTED_POTENTIAL:
1167      case esm_TAYLOR_SHIFTED:
1168 +    case esm_EWALD_FULL:
1169        if (i_is_Charge)
1170          self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge));      
1171        if (i_is_Dipole)
# Line 1192 | Line 1193 | namespace OpenMD {
1193    }
1194  
1195  
1196 <  void Electrostatic::ReciprocalSpaceSum () {
1196 >  void Electrostatic::ReciprocalSpaceSum(RealType& pot) {
1197      
1198      RealType kPot = 0.0;
1199      RealType kVir = 0.0;
# Line 1217 | Line 1218 | namespace OpenMD {
1218      Vector3d box = hmat.diagonals();
1219      RealType boxMax = box.max();
1220      
1221 <    //int kMax = int(pow(dampingAlpha_,2)*cutoffRadius_ * boxMax / M_PI);
1222 <    const int kMax = 5;
1221 >    //int kMax = int(2.0 * M_PI / (pow(dampingAlpha_,2)*cutoffRadius_ * boxMax) );
1222 >    int kMax = 7;
1223      int kSqMax = kMax*kMax + 2;
1224      
1225      int kLimit = kMax+1;
# Line 1238 | Line 1239 | namespace OpenMD {
1239      
1240      // Calculate and store exponential factors  
1241      
1242 <    vector<vector<Vector3d> > eCos;
1243 <    vector<vector<Vector3d> > eSin;
1242 >    vector<vector<RealType> > elc;
1243 >    vector<vector<RealType> > emc;
1244 >    vector<vector<RealType> > enc;
1245 >    vector<vector<RealType> > els;
1246 >    vector<vector<RealType> > ems;
1247 >    vector<vector<RealType> > ens;
1248 >
1249      
1250      int nMax = info_->getNAtoms();
1251      
1252 <    eCos.resize(kLimit+1);
1253 <    eSin.resize(kLimit+1);
1252 >    elc.resize(kLimit+1);
1253 >    emc.resize(kLimit+1);
1254 >    enc.resize(kLimit+1);
1255 >    els.resize(kLimit+1);
1256 >    ems.resize(kLimit+1);
1257 >    ens.resize(kLimit+1);
1258 >
1259      for (int j = 0; j < kLimit+1; j++) {
1260 <      eCos[j].resize(nMax);
1261 <      eSin[j].resize(nMax);
1260 >      elc[j].resize(nMax);
1261 >      emc[j].resize(nMax);
1262 >      enc[j].resize(nMax);
1263 >      els[j].resize(nMax);
1264 >      ems[j].resize(nMax);
1265 >      ens[j].resize(nMax);
1266      }
1267      
1268      Vector3d t( 2.0 * M_PI );
1269      t.Vdiv(t, box);
1270 +
1271      
1272      SimInfo::MoleculeIterator mi;
1273      Molecule::AtomIterator ai;
1274      int i;
1275      Vector3d r;
1276      Vector3d tt;
1261 –    Vector3d w;
1262 –    Vector3d u;
1277      
1278      for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1279           mol = info_->nextMolecule(mi)) {
# Line 1271 | Line 1285 | namespace OpenMD {
1285          info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r);
1286          
1287          tt.Vmul(t, r);
1288 +
1289 +        elc[1][i] = 1.0;
1290 +        emc[1][i] = 1.0;
1291 +        enc[1][i] = 1.0;
1292 +        els[1][i] = 0.0;
1293 +        ems[1][i] = 0.0;
1294 +        ens[1][i] = 0.0;
1295 +
1296 +        elc[2][i] = cos(tt.x());
1297 +        emc[2][i] = cos(tt.y());
1298 +        enc[2][i] = cos(tt.z());
1299 +        els[2][i] = sin(tt.x());
1300 +        ems[2][i] = sin(tt.y());
1301 +        ens[2][i] = sin(tt.z());
1302          
1275 –        eCos[1][i] = Vector3d(1.0, 1.0, 1.0);
1276 –        eSin[1][i] = Vector3d(0.0, 0.0, 0.0);
1277 –        eCos[2][i] = Vector3d(cos(tt.x()), cos(tt.y()), cos(tt.z()));
1278 –        eSin[2][i] = Vector3d(sin(tt.x()), sin(tt.y()), sin(tt.z()));
1279 –        u = 2.0 * eCos[1][i];
1280 –        eCos[3][i].Vmul(u, eCos[2][i]);
1281 –        eSin[3][i].Vmul(u, eSin[2][i]);      
1282 –        
1303          for(int l = 3; l <= kLimit; l++) {
1304 <          w.Vmul(u, eCos[l-1][i]);
1305 <          eCos[l][i] = w - eCos[l-2][i];
1306 <          w.Vmul(u, eSin[l-1][i]);
1307 <          eSin[l][i] = w - eSin[l-2][i];
1304 >          elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i];
1305 >          emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i];
1306 >          enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i];
1307 >          els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i];
1308 >          ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i];
1309 >          ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i];
1310          }
1311        }
1312      }
# Line 1328 | Line 1350 | namespace OpenMD {
1350      std::vector<RealType> qks(nMax, 0.0);
1351      std::vector<Vector3d> dxk(nMax, V3Zero);
1352      std::vector<Vector3d> qxk(nMax, V3Zero);
1353 <    
1353 >    RealType rl, rm, rn;
1354 >    Vector3d kVec;
1355 >    Vector3d Qk;
1356 >    Mat3x3d k2;
1357 >    RealType ckcs, ckss, dkcs, dkss, qkcs, qkss;
1358 >    int atid;
1359 >    ElectrostaticAtomData data;
1360 >    RealType C, dk, qk;
1361 >    Vector3d D;
1362 >    Mat3x3d  Q;
1363 >
1364      int mMin = kLimit;
1365      int nMin = kLimit + 1;
1366      for (int l = 1; l <= kLimit; l++) {
1367 <      int ll =l - 1;
1368 <      RealType rl = xcl * float(ll);
1367 >      int ll = l - 1;
1368 >      rl = xcl * float(ll);
1369        for (int mmm = mMin; mmm <= kLim2; mmm++) {
1370          int mm = mmm - kLimit;
1371          int m = abs(mm) + 1;
1372 <        RealType rm = ycl * float(mm);
1372 >        rm = ycl * float(mm);
1373          // Set temporary products of exponential terms
1374          for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1375               mol = info_->nextMolecule(mi)) {
# Line 1346 | Line 1378 | namespace OpenMD {
1378              
1379              i = atom->getLocalIndex();
1380              if(mm < 0) {
1381 <              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1382 <                     + eSin[l][i].x()*eSin[m][i].y();
1351 <              slm[i] = eCos[l][i].x()*eCos[m][i].y()
1352 <                     - eSin[m][i].y()*eCos[l][i].x();
1381 >              clm[i]=elc[l][i]*emc[m][i]+els[l][i]*ems[m][i];
1382 >              slm[i]=els[l][i]*emc[m][i]-ems[m][i]*elc[l][i];
1383              } else {
1384 <              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1385 <                     - eSin[l][i].x()*eSin[m][i].y();
1356 <              slm[i] = eSin[l][i].x()*eCos[m][i].y()
1357 <                     + eSin[m][i].y()*eCos[l][i].x();
1384 >              clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i];
1385 >              slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i];
1386              }
1387            }
1388          }
1389          for (int nnn = nMin; nnn <= kLim2; nnn++) {
1390            int nn = nnn - kLimit;          
1391            int n = abs(nn) + 1;
1392 <          RealType rn = zcl * float(nn);
1392 >          rn = zcl * float(nn);
1393            // Test on magnitude of k vector:
1394            int kk=ll*ll + mm*mm + nn*nn;
1395            if(kk <= kSqLim) {
1396 <            Vector3d kVec = Vector3d(rl, rm, rn);
1397 <            Mat3x3d  k2 = outProduct(kVec, kVec);
1396 >            kVec = Vector3d(rl, rm, rn);
1397 >            k2 = outProduct(kVec, kVec);
1398              // Calculate exp(ikr) terms
1399              for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1400                   mol = info_->nextMolecule(mi)) {
# Line 1375 | Line 1403 | namespace OpenMD {
1403                  i = atom->getLocalIndex();
1404                  
1405                  if (nn < 0) {
1406 <                  ckr[i]=clm[i]*eCos[n][i].z()+slm[i]*eSin[n][i].z();
1407 <                  skr[i]=slm[i]*eCos[n][i].z()-clm[i]*eSin[n][i].z();
1406 >                  ckr[i]=clm[i]*enc[n][i]+slm[i]*ens[n][i];
1407 >                  skr[i]=slm[i]*enc[n][i]-clm[i]*ens[n][i];
1408 >
1409                  } else {
1410 <                  ckr[i]=clm[i]*eCos[n][i].z()-slm[i]*eSin[n][i].z();
1411 <                  skr[i]=slm[i]*eCos[n][i].z()+clm[i]*eSin[n][i].z();
1410 >                  ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i];
1411 >                  skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i];
1412                  }
1413                }
1414              }
# Line 1390 | Line 1419 | namespace OpenMD {
1419                   mol = info_->nextMolecule(mi)) {
1420                for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1421                    atom = mol->nextAtom(ai)) {
1422 <                i = atom->getGlobalIndex();
1422 >                i = atom->getLocalIndex();
1423                  int atid = atom->getAtomType()->getIdent();
1424 <                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1424 >                data = ElectrostaticMap[Etids[atid]];
1425                                
1426                  if (data.is_Charge) {
1427 <                  RealType C = data.fixedCharge;
1427 >                  C = data.fixedCharge;
1428                    if (atom->isFluctuatingCharge()) C += atom->getFlucQPos();
1429                    ckc[i] = C * ckr[i];
1430 <                  cks[i] = C * cks[i];
1430 >                  cks[i] = C * skr[i];
1431                  }
1432                  
1433                  if (data.is_Dipole) {
1434 <                  Vector3d D = atom->getDipole() * mPoleConverter;
1435 <                  RealType dk = dot(kVec, D);
1436 <                  dxk[i] = cross(kVec, D);
1434 >                  D = atom->getDipole() * mPoleConverter;
1435 >                  dk = dot(D, kVec);
1436 >                  dxk[i] = cross(D, kVec);
1437                    dkc[i] = dk * ckr[i];
1438                    dks[i] = dk * skr[i];
1439                  }
1440                  if (data.is_Quadrupole) {
1441 <                  Mat3x3d Q = atom->getQuadrupole();
1442 <                  Q *= mPoleConverter;
1443 <                  RealType qk = -( Q * k2 ).trace();
1444 <                  qxk[i] = -2.0 * cross(k2, Q);
1441 >                  Q = atom->getQuadrupole() * mPoleConverter;
1442 >                  Qk = Q * kVec;                  
1443 >                  qk = dot(kVec, Qk);
1444 >                  qxk[i] = -cross(kVec, Qk);
1445                    qkc[i] = qk * ckr[i];
1446                    qks[i] = qk * skr[i];
1447                  }              
1448                }
1449              }
1450 <            
1450 >
1451              // calculate vector sums
1452              
1453 <            RealType ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1454 <            RealType ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1455 <            RealType dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1456 <            RealType dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1457 <            RealType qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1458 <            RealType qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1453 >            ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1454 >            ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1455 >            dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1456 >            dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1457 >            qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1458 >            qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1459              
1460   #ifdef IS_MPI
1461              MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE,
# Line 1444 | Line 1473 | namespace OpenMD {
1473   #endif        
1474              
1475              // Accumulate potential energy and virial contribution:
1476 +
1477 +            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1478 +                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs));
1479 +
1480 +            kVir += 2.0 * rvol  * AK[kk]*(ckcs*ckcs+ckss*ckss
1481 +                                          +4.0*(ckss*dkcs-ckcs*dkss)
1482 +                                          +3.0*(dkcs*dkcs+dkss*dkss)
1483 +                                          -6.0*(ckss*qkss+ckcs*qkcs)
1484 +                                          +8.0*(dkss*qkcs-dkcs*qkss)
1485 +                                          +5.0*(qkss*qkss+qkcs*qkcs));
1486              
1448 –            //cerr << "l, m, n = " << l << " " << m << " " << n << "\n";
1449 –            cerr << "kVec = " << kVec << "\n";
1450 –            cerr << "ckss = " << ckss << " ckcs = " << ckcs << "\n";
1451 –            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1452 –                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkss));
1453 –            //cerr << "kspace pot = " << kPot << "\n";
1454 –            kVir -= 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss
1455 –                                         +4.0*(ckss*dkcs-ckcs*dkss)
1456 –                                         +3.0*(dkcs*dkcs+dkss*dkss)
1457 –                                         -6.0*(ckss*qkss+ckcs*qkcs)
1458 –                                         +8.0*(dkss*qkcs-dkcs*qkss)
1459 –                                         +5.0*(qkss*qkss+qkcs*qkcs));
1460 –            
1487              // Calculate force and torque for each site:
1488              
1489              for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
# Line 1466 | Line 1492 | namespace OpenMD {
1492                    atom = mol->nextAtom(ai)) {
1493                  
1494                  i = atom->getLocalIndex();
1495 <                int atid = atom->getAtomType()->getIdent();
1496 <                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1497 <                
1495 >                atid = atom->getAtomType()->getIdent();
1496 >                data = ElectrostaticMap[Etids[atid]];
1497 >
1498                  RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs)
1499 <                                        - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss));
1499 >                                     - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss));
1500                  RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
1501                                           -ckr[i]*(ckss+dkcs-qkss));
1502 <                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)+
1503 <                                             skr[i]*(ckss+dkcs-qkss));
1504 <                
1479 <                
1502 >                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
1503 >                                            +skr[i]*(ckss+dkcs-qkss));
1504 >              
1505                  atom->addFrc( 4.0 * rvol * qfrc * kVec );
1506                  
1507                  if (data.is_Dipole) {
# Line 1489 | Line 1514 | namespace OpenMD {
1514              }
1515            }
1516          }
1517 +        nMin = 1;
1518        }
1519 +      mMin = 1;
1520      }
1521 +    pot += kPot;  
1522    }
1523   }

Diff Legend

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