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 1921 by gezelter, Thu Aug 1 18:23:07 2013 UTC vs.
Revision 1925 by gezelter, Wed Aug 7 15:24:16 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 289 | Line 291 | namespace OpenMD {
291      db0c_4 =          3.0*b2c  - 6.0*r2*b3c     + r2*r2*b4c;
292      db0c_5 =                    -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c;  
293  
294 <    if (summationMethod_ == esm_EWALD_FULL) {
293 <      selfMult1_ *= 2.0;
294 <      selfMult2_ *= 2.0;
295 <      selfMult4_ *= 2.0;
296 <    } else {
294 >    if (summationMethod_ != esm_EWALD_FULL) {
295        selfMult1_ -= b0c;
296        selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) /  3.0;
297        selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0;
# Line 1198 | Line 1196 | namespace OpenMD {
1196    }
1197  
1198  
1199 <  void Electrostatic::ReciprocalSpaceSum(potVec& pot) {
1199 >  void Electrostatic::ReciprocalSpaceSum(RealType& pot) {
1200      
1201      RealType kPot = 0.0;
1202      RealType kVir = 0.0;
# Line 1223 | Line 1221 | namespace OpenMD {
1221      Vector3d box = hmat.diagonals();
1222      RealType boxMax = box.max();
1223      
1226 –    cerr << "da = " << dampingAlpha_ << " rc = " << cutoffRadius_ << "\n";
1227 –    cerr << "boxMax = " << boxMax << "\n";
1224      //int kMax = int(2.0 * M_PI / (pow(dampingAlpha_,2)*cutoffRadius_ * boxMax) );
1225 <    int kMax = 5;
1230 <    cerr << "kMax = " << kMax << "\n";
1225 >    int kMax = 7;
1226      int kSqMax = kMax*kMax + 2;
1227      
1228      int kLimit = kMax+1;
# Line 1247 | Line 1242 | namespace OpenMD {
1242      
1243      // Calculate and store exponential factors  
1244      
1245 <    vector<vector<Vector3d> > eCos;
1246 <    vector<vector<Vector3d> > eSin;
1245 >    vector<vector<RealType> > elc;
1246 >    vector<vector<RealType> > emc;
1247 >    vector<vector<RealType> > enc;
1248 >    vector<vector<RealType> > els;
1249 >    vector<vector<RealType> > ems;
1250 >    vector<vector<RealType> > ens;
1251 >
1252      
1253      int nMax = info_->getNAtoms();
1254      
1255 <    eCos.resize(kLimit+1);
1256 <    eSin.resize(kLimit+1);
1255 >    elc.resize(kLimit+1);
1256 >    emc.resize(kLimit+1);
1257 >    enc.resize(kLimit+1);
1258 >    els.resize(kLimit+1);
1259 >    ems.resize(kLimit+1);
1260 >    ens.resize(kLimit+1);
1261 >
1262      for (int j = 0; j < kLimit+1; j++) {
1263 <      eCos[j].resize(nMax);
1264 <      eSin[j].resize(nMax);
1263 >      elc[j].resize(nMax);
1264 >      emc[j].resize(nMax);
1265 >      enc[j].resize(nMax);
1266 >      els[j].resize(nMax);
1267 >      ems[j].resize(nMax);
1268 >      ens[j].resize(nMax);
1269      }
1270      
1271      Vector3d t( 2.0 * M_PI );
1272      t.Vdiv(t, box);
1273 +
1274      
1275      SimInfo::MoleculeIterator mi;
1276      Molecule::AtomIterator ai;
1277      int i;
1278      Vector3d r;
1279      Vector3d tt;
1270 –    Vector3d w;
1271 –    Vector3d u;
1272 –    Vector3d a;
1273 –    Vector3d b;
1280      
1281      for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1282           mol = info_->nextMolecule(mi)) {
# Line 1281 | Line 1287 | namespace OpenMD {
1287          r = atom->getPos();
1288          info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r);
1289          
1284 –        // Shift so that all coordinates are in the range [0,L]:
1285 –
1286 –        r += box/2.0;
1287 –
1290          tt.Vmul(t, r);
1291  
1292 <        //cerr << "tt = " << tt << "\n";
1292 >        elc[1][i] = 1.0;
1293 >        emc[1][i] = 1.0;
1294 >        enc[1][i] = 1.0;
1295 >        els[1][i] = 0.0;
1296 >        ems[1][i] = 0.0;
1297 >        ens[1][i] = 0.0;
1298 >
1299 >        elc[2][i] = cos(tt.x());
1300 >        emc[2][i] = cos(tt.y());
1301 >        enc[2][i] = cos(tt.z());
1302 >        els[2][i] = sin(tt.x());
1303 >        ems[2][i] = sin(tt.y());
1304 >        ens[2][i] = sin(tt.z());
1305          
1292 –        eCos[1][i] = Vector3d(1.0, 1.0, 1.0);
1293 –        eSin[1][i] = Vector3d(0.0, 0.0, 0.0);
1294 –        eCos[2][i] = Vector3d(cos(tt.x()), cos(tt.y()), cos(tt.z()));
1295 –        eSin[2][i] = Vector3d(sin(tt.x()), sin(tt.y()), sin(tt.z()));
1296 –        u = eCos[2][i];
1297 –        w = eSin[2][i];
1298 –        
1306          for(int l = 3; l <= kLimit; l++) {
1307 <          a.Vmul(eCos[l-1][i], u);
1308 <          b.Vmul(eSin[l-1][i], w);
1309 <          eCos[l][i] = a - b;
1310 <          a.Vmul(eSin[l-1][i], u);
1311 <          b.Vmul(eCos[l-1][i], w);
1312 <          eSin[l][i] = a + b;
1307 >          elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i];
1308 >          emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i];
1309 >          enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i];
1310 >          els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i];
1311 >          ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i];
1312 >          ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i];
1313          }
1314        }
1315      }
# Line 1346 | Line 1353 | namespace OpenMD {
1353      std::vector<RealType> qks(nMax, 0.0);
1354      std::vector<Vector3d> dxk(nMax, V3Zero);
1355      std::vector<Vector3d> qxk(nMax, V3Zero);
1356 <    
1356 >    RealType rl, rm, rn;
1357 >    Vector3d kVec;
1358 >    Vector3d Qk;
1359 >    Mat3x3d k2;
1360 >    RealType ckcs, ckss, dkcs, dkss, qkcs, qkss;
1361 >    int atid;
1362 >    ElectrostaticAtomData data;
1363 >    RealType C, dk, qk;
1364 >    Vector3d D;
1365 >    Mat3x3d  Q;
1366 >
1367      int mMin = kLimit;
1368      int nMin = kLimit + 1;
1369      for (int l = 1; l <= kLimit; l++) {
1370 <      int ll =l - 1;
1371 <      RealType rl = xcl * float(ll);
1370 >      int ll = l - 1;
1371 >      rl = xcl * float(ll);
1372        for (int mmm = mMin; mmm <= kLim2; mmm++) {
1373          int mm = mmm - kLimit;
1374          int m = abs(mm) + 1;
1375 <        RealType rm = ycl * float(mm);
1375 >        rm = ycl * float(mm);
1376          // Set temporary products of exponential terms
1377          for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1378               mol = info_->nextMolecule(mi)) {
# Line 1364 | Line 1381 | namespace OpenMD {
1381              
1382              i = atom->getLocalIndex();
1383              if(mm < 0) {
1384 <              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1385 <                     + 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();
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              } else {
1387 <              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1388 <                     - 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();
1387 >              clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i];
1388 >              slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i];
1389              }
1390            }
1391          }
1392          for (int nnn = nMin; nnn <= kLim2; nnn++) {
1393            int nn = nnn - kLimit;          
1394            int n = abs(nn) + 1;
1395 <          RealType rn = zcl * float(nn);
1395 >          rn = zcl * float(nn);
1396            // Test on magnitude of k vector:
1397            int kk=ll*ll + mm*mm + nn*nn;
1398            if(kk <= kSqLim) {
1399 <            Vector3d kVec = Vector3d(rl, rm, rn);
1400 <            Mat3x3d  k2 = outProduct(kVec, kVec);
1399 >            kVec = Vector3d(rl, rm, rn);
1400 >            k2 = outProduct(kVec, kVec);
1401              // Calculate exp(ikr) terms
1402              for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1403                   mol = info_->nextMolecule(mi)) {
# Line 1393 | Line 1406 | namespace OpenMD {
1406                  i = atom->getLocalIndex();
1407                  
1408                  if (nn < 0) {
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                  } else {
1413 <                  ckr[i]=clm[i]*eCos[n][i].z()-slm[i]*eSin[n][i].z();
1414 <                  skr[i]=slm[i]*eCos[n][i].z()+clm[i]*eSin[n][i].z();
1413 >                  ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i];
1414 >                  skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i];
1415                  }
1416                }
1417              }
# Line 1410 | Line 1424 | namespace OpenMD {
1424                    atom = mol->nextAtom(ai)) {
1425                  i = atom->getLocalIndex();
1426                  int atid = atom->getAtomType()->getIdent();
1427 <                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1427 >                data = ElectrostaticMap[Etids[atid]];
1428                                
1429                  if (data.is_Charge) {
1430 <                  RealType C = data.fixedCharge;
1430 >                  C = data.fixedCharge;
1431                    if (atom->isFluctuatingCharge()) C += atom->getFlucQPos();
1432                    ckc[i] = C * ckr[i];
1433                    cks[i] = C * skr[i];
1434                  }
1435                  
1436                  if (data.is_Dipole) {
1437 <                  Vector3d D = atom->getDipole() * mPoleConverter;
1438 <                  RealType dk = dot(kVec, D);
1439 <                  dxk[i] = cross(kVec, D);
1437 >                  D = atom->getDipole() * mPoleConverter;
1438 >                  dk = dot(D, kVec);
1439 >                  dxk[i] = cross(D, kVec);
1440                    dkc[i] = dk * ckr[i];
1441                    dks[i] = dk * skr[i];
1442                  }
1443                  if (data.is_Quadrupole) {
1444 <                  Mat3x3d Q = atom->getQuadrupole();
1444 >                  Q = atom->getQuadrupole();
1445                    Q *= mPoleConverter;
1446 <                  RealType qk = -( Q * k2 ).trace();
1447 <                  qxk[i] = -2.0 * cross(k2, Q);
1446 >                  Qk = Q * kVec;
1447 >                  qk = dot(kVec, Qk);
1448 >                  qxk[i] = cross(kVec, Qk);
1449                    qkc[i] = qk * ckr[i];
1450                    qks[i] = qk * skr[i];
1451                  }              
# Line 1439 | Line 1454 | namespace OpenMD {
1454  
1455              // calculate vector sums
1456              
1457 <            RealType ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1458 <            RealType ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1459 <            RealType dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1460 <            RealType dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1461 <            RealType qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1462 <            RealType qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1448 <
1457 >            ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1458 >            ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1459 >            dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1460 >            dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1461 >            qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1462 >            qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1463              
1464   #ifdef IS_MPI
1465              MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE,
# Line 1462 | Line 1476 | namespace OpenMD {
1476                                        MPI::SUM);
1477   #endif        
1478              
1479 <            // Accumulate potential energy and virial contribution:          
1479 >            // Accumulate potential energy and virial contribution:
1480  
1481 <            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1482 <                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkss));
1481 >            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1482 >                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs));
1483  
1484 <            kVir -= 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss
1485 <                                         +4.0*(ckss*dkcs-ckcs*dkss)
1486 <                                         +3.0*(dkcs*dkcs+dkss*dkss)
1487 <                                         -6.0*(ckss*qkss+ckcs*qkcs)
1488 <                                         +8.0*(dkss*qkcs-dkcs*qkss)
1489 <                                         +5.0*(qkss*qkss+qkcs*qkcs));
1484 >            kVir += 2.0 * rvol  * AK[kk]*(ckcs*ckcs+ckss*ckss
1485 >                                          +4.0*(ckss*dkcs-ckcs*dkss)
1486 >                                          +3.0*(dkcs*dkcs+dkss*dkss)
1487 >                                          -6.0*(ckss*qkss+ckcs*qkcs)
1488 >                                          +8.0*(dkss*qkcs-dkcs*qkss)
1489 >                                          +5.0*(qkss*qkss+qkcs*qkcs));
1490              
1491              // Calculate force and torque for each site:
1492              
# Line 1482 | Line 1496 | namespace OpenMD {
1496                    atom = mol->nextAtom(ai)) {
1497                  
1498                  i = atom->getLocalIndex();
1499 <                int atid = atom->getAtomType()->getIdent();
1500 <                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1501 <                
1499 >                atid = atom->getAtomType()->getIdent();
1500 >                data = ElectrostaticMap[Etids[atid]];
1501 >
1502                  RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs)
1503 <                                        - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss));
1503 >                                     - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss));
1504                  RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
1505                                           -ckr[i]*(ckss+dkcs-qkss));
1506 <                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)+
1507 <                                             skr[i]*(ckss+dkcs-qkss));
1508 <                
1495 <                
1506 >                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
1507 >                                             +skr[i]*(ckss+dkcs-qkss));
1508 >              
1509                  atom->addFrc( 4.0 * rvol * qfrc * kVec );
1510                  
1511                  if (data.is_Dipole) {
# Line 1505 | Line 1518 | namespace OpenMD {
1518              }
1519            }
1520          }
1521 +        nMin = 1;
1522        }
1523 +      mMin = 1;
1524      }
1525 <    cerr << "kPot = " << kPot << "\n";
1511 <    pot[ELECTROSTATIC_FAMILY] += kPot;  
1525 >    pot += kPot;  
1526    }
1527   }

Diff Legend

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