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 1900 by gezelter, Fri Jul 12 17:38:06 2013 UTC vs.
Revision 1915 by gezelter, Mon Jul 29 17:55:17 2013 UTC

# Line 44 | Line 44
44   #include <string.h>
45  
46   #include <cmath>
47 + #include <numeric>
48   #include "nonbonded/Electrostatic.hpp"
49   #include "utils/simError.h"
50   #include "types/NonBondedInteractionType.hpp"
# Line 55 | Line 56
56   #include "utils/PhysicalConstants.hpp"
57   #include "math/erfc.hpp"
58   #include "math/SquareMatrix.hpp"
59 + #include "primitives/Molecule.hpp"
60  
61 +
62   namespace OpenMD {
63    
64    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
# Line 191 | Line 194 | namespace OpenMD {
194        simError();
195      }
196            
197 <    if (screeningMethod_ == DAMPED) {      
197 >    if (screeningMethod_ == DAMPED || summationMethod_ == esm_EWALD_FULL) {
198        if (!simParams_->haveDampingAlpha()) {
199          // first set a cutoff dependent alpha value
200          // we assume alpha depends linearly with rcut from 0 to 20.5 ang
201          dampingAlpha_ = 0.425 - cutoffRadius_* 0.02;
202 <        if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;
200 <        
202 >        if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;        
203          // throw warning
204          sprintf( painCave.errMsg,
205                   "Electrostatic::initialize: dampingAlpha was not specified in the\n"
# Line 213 | Line 215 | namespace OpenMD {
215        haveDampingAlpha_ = true;
216      }
217  
218 +
219      Etypes.clear();
220      Etids.clear();
221      FQtypes.clear();
# Line 262 | Line 265 | namespace OpenMD {
265        b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2;
266        b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2;
267        b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2;
265 –      //selfMult1_ = - 2.0 * a2 * invArootPi;
266 –      //selfMult2_ = - 4.0 * a2 * a2 * invArootPi / 3.0;
267 –      //selfMult4_ = - 8.0 * a2 * a2 * a2 * invArootPi / 5.0;
268        // Half the Smith self piece:
269        selfMult1_ = - a2 * invArootPi;
270        selfMult2_ = - 2.0 * a2 * a2 * invArootPi / 3.0;
# Line 288 | Line 288 | namespace OpenMD {
288      db0c_3 =          3.0*r*b2c  - r2*r*b3c;
289      db0c_4 =          3.0*b2c  - 6.0*r2*b3c     + r2*r2*b4c;
290      db0c_5 =                    -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c;  
291 –    
292 –    selfMult1_ -= b0c;
293 –    selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) /  3.0;
294 –    selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0;
291  
292 +    if (summationMethod_ != esm_EWALD_FULL) {
293 +      selfMult1_ -= b0c;
294 +      selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) /  3.0;
295 +      selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0;
296 +    }
297 +
298      // working variables for the splines:
299      RealType ri, ri2;
300      RealType b0, b1, b2, b3, b4, b5;
# Line 328 | Line 330 | namespace OpenMD {
330      vector<RealType> v31v, v32v;
331      vector<RealType> v41v, v42v, v43v;
332  
331 –    /*
332 –    vector<RealType> dv01v;
333 –    vector<RealType> dv11v;
334 –    vector<RealType> dv21v, dv22v;
335 –    vector<RealType> dv31v, dv32v;
336 –    vector<RealType> dv41v, dv42v, dv43v;
337 –    */
338 –
333      for (int i = 1; i < np_ + 1; i++) {
334        r = RealType(i) * dx;
335        rv.push_back(r);
# Line 499 | Line 493 | namespace OpenMD {
493  
494        case esm_SWITCHING_FUNCTION:
495        case esm_HARD:
496 +      case esm_EWALD_FULL:
497  
498          v01 = f;
499          v11 = g;
# Line 557 | Line 552 | namespace OpenMD {
552  
553          break;
554                  
560 –      case esm_EWALD_FULL:
555        case esm_EWALD_PME:
556        case esm_EWALD_SPME:
557        default :
# Line 586 | Line 580 | namespace OpenMD {
580        v41v.push_back(v41);
581        v42v.push_back(v42);
582        v43v.push_back(v43);
589 –      /*
590 –      dv01v.push_back(dv01);
591 –      dv11v.push_back(dv11);
592 –      dv21v.push_back(dv21);
593 –      dv22v.push_back(dv22);
594 –      dv31v.push_back(dv31);
595 –      dv32v.push_back(dv32);      
596 –      dv41v.push_back(dv41);
597 –      dv42v.push_back(dv42);
598 –      dv43v.push_back(dv43);
599 –      */
583      }
584  
585      // construct the spline structures and fill them with the values we've
# Line 621 | Line 604 | namespace OpenMD {
604      v43s = new CubicSpline();
605      v43s->addPoints(rv, v43v);
606  
624 –    /*
625 –    dv01s = new CubicSpline();
626 –    dv01s->addPoints(rv, dv01v);
627 –    dv11s = new CubicSpline();
628 –    dv11s->addPoints(rv, dv11v);
629 –    dv21s = new CubicSpline();
630 –    dv21s->addPoints(rv, dv21v);
631 –    dv22s = new CubicSpline();
632 –    dv22s->addPoints(rv, dv22v);
633 –    dv31s = new CubicSpline();
634 –    dv31s->addPoints(rv, dv31v);
635 –    dv32s = new CubicSpline();
636 –    dv32s->addPoints(rv, dv32v);
637 –    dv41s = new CubicSpline();
638 –    dv41s->addPoints(rv, dv41v);
639 –    dv42s = new CubicSpline();
640 –    dv42s->addPoints(rv, dv42v);
641 –    dv43s = new CubicSpline();
642 –    dv43s->addPoints(rv, dv43v);
643 –    */
644 –
607      haveElectroSplines_ = true;
608  
609      initialized_ = true;
# Line 1202 | Line 1164 | namespace OpenMD {
1164        
1165      case esm_SHIFTED_FORCE:
1166      case esm_SHIFTED_POTENTIAL:
1167 +    case esm_TAYLOR_SHIFTED:
1168        if (i_is_Charge)
1169          self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge));      
1170        if (i_is_Dipole)
# Line 1227 | Line 1190 | namespace OpenMD {
1190      // cases.
1191      return 12.0;
1192    }
1193 +
1194 +
1195 +  void Electrostatic::ReciprocalSpaceSum () {
1196 +    
1197 +    RealType kPot = 0.0;
1198 +    RealType kVir = 0.0;
1199 +    
1200 +    const RealType mPoleConverter = 0.20819434; // converts from the
1201 +                                                // internal units of
1202 +                                                // Debye (for dipoles)
1203 +                                                // or Debye-angstroms
1204 +                                                // (for quadrupoles) to
1205 +                                                // electron angstroms or
1206 +                                                // electron-angstroms^2
1207 +    
1208 +    const RealType eConverter = 332.0637778; // convert the
1209 +                                             // Charge-Charge
1210 +                                             // electrostatic
1211 +                                             // interactions into kcal /
1212 +                                             // mol assuming distances
1213 +                                             // are measured in
1214 +                                             // angstroms.
1215 +
1216 +    Mat3x3d hmat = info_->getSnapshotManager()->getCurrentSnapshot()->getHmat();
1217 +    Vector3d box = hmat.diagonals();
1218 +    RealType boxMax = box.max();
1219 +    
1220 +    //int kMax = int(pow(dampingAlpha_,2)*cutoffRadius_ * boxMax / M_PI);
1221 +    const int kMax = 5;
1222 +    int kSqMax = kMax*kMax + 2;
1223 +    
1224 +    int kLimit = kMax+1;
1225 +    int kLim2 = 2*kMax+1;
1226 +    int kSqLim = kSqMax;
1227 +    
1228 +    vector<RealType> AK(kSqLim+1, 0.0);
1229 +    RealType xcl = 2.0 * M_PI / box.x();
1230 +    RealType ycl = 2.0 * M_PI / box.y();
1231 +    RealType zcl = 2.0 * M_PI / box.z();
1232 +    RealType rcl = 2.0 * M_PI / boxMax;
1233 +    RealType rvol = 2.0 * M_PI /(box.x() * box.y() * box.z());
1234 +    
1235 +    if(dampingAlpha_ < 1.0e-12) return;
1236 +    
1237 +    RealType ralph = -0.25/pow(dampingAlpha_,2);
1238 +    
1239 +    // Calculate and store exponential factors  
1240 +    
1241 +    vector<vector<Vector3d> > eCos;
1242 +    vector<vector<Vector3d> > eSin;
1243 +    
1244 +    int nMax = info_->getNAtoms();
1245 +    
1246 +    eCos.resize(kLimit+1);
1247 +    eSin.resize(kLimit+1);
1248 +    for (int j = 0; j < kLimit+1; j++) {
1249 +      eCos[j].resize(nMax);
1250 +      eSin[j].resize(nMax);
1251 +    }
1252 +    
1253 +    Vector3d t( 2.0 * M_PI );
1254 +    t.Vdiv(t, box);
1255 +    
1256 +    SimInfo::MoleculeIterator mi;
1257 +    Molecule::AtomIterator ai;
1258 +    int i;
1259 +    Vector3d r;
1260 +    Vector3d tt;
1261 +    Vector3d w;
1262 +    Vector3d u;
1263 +    
1264 +    for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1265 +         mol = info_->nextMolecule(mi)) {
1266 +      for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1267 +          atom = mol->nextAtom(ai)) {  
1268 +        
1269 +        i = atom->getLocalIndex();
1270 +        r = atom->getPos();
1271 +        info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r);
1272 +        
1273 +        tt.Vmul(t, r);
1274 +        
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 +        
1283 +        for(int l = 3; l <= kLimit; l++) {
1284 +          w.Vmul(u, eCos[l-1][i]);
1285 +          eCos[l][i] = w - eCos[l-2][i];
1286 +          w.Vmul(u, eSin[l-1][i]);
1287 +          eSin[l][i] = w - eSin[l-2][i];
1288 +        }
1289 +      }
1290 +    }
1291 +    
1292 +    // Calculate and store AK coefficients:
1293 +    
1294 +    RealType eksq = 1.0;
1295 +    RealType expf = 0.0;
1296 +    if (ralph < 0.0) expf = exp(ralph*rcl*rcl);
1297 +    for (i = 1; i <= kSqLim; i++) {
1298 +      RealType rksq = float(i)*rcl*rcl;
1299 +      eksq = expf*eksq;
1300 +      AK[i] = eConverter * eksq/rksq;
1301 +    }
1302 +    
1303 +    /*
1304 +     * Loop over all k vectors k = 2 pi (ll/Lx, mm/Ly, nn/Lz)
1305 +     * the values of ll, mm and nn are selected so that the symmetry of
1306 +     * reciprocal lattice is taken into account i.e. the following
1307 +     * rules apply.
1308 +     *
1309 +     * ll ranges over the values 0 to kMax only.
1310 +     *
1311 +     * mm ranges over 0 to kMax when ll=0 and over
1312 +     *            -kMax to kMax otherwise.
1313 +     * nn ranges over 1 to kMax when ll=mm=0 and over
1314 +     *            -kMax to kMax otherwise.
1315 +     *
1316 +     * Hence the result of the summation must be doubled at the end.    
1317 +     */
1318 +    
1319 +    std::vector<RealType> clm(nMax, 0.0);
1320 +    std::vector<RealType> slm(nMax, 0.0);
1321 +    std::vector<RealType> ckr(nMax, 0.0);
1322 +    std::vector<RealType> skr(nMax, 0.0);
1323 +    std::vector<RealType> ckc(nMax, 0.0);
1324 +    std::vector<RealType> cks(nMax, 0.0);
1325 +    std::vector<RealType> dkc(nMax, 0.0);
1326 +    std::vector<RealType> dks(nMax, 0.0);
1327 +    std::vector<RealType> qkc(nMax, 0.0);
1328 +    std::vector<RealType> qks(nMax, 0.0);
1329 +    std::vector<Vector3d> dxk(nMax, V3Zero);
1330 +    std::vector<Vector3d> qxk(nMax, V3Zero);
1331 +    
1332 +    int mMin = kLimit;
1333 +    int nMin = kLimit + 1;
1334 +    for (int l = 1; l <= kLimit; l++) {
1335 +      int ll =l - 1;
1336 +      RealType rl = xcl * float(ll);
1337 +      for (int mmm = mMin; mmm <= kLim2; mmm++) {
1338 +        int mm = mmm - kLimit;
1339 +        int m = abs(mm) + 1;
1340 +        RealType rm = ycl * float(mm);
1341 +        // Set temporary products of exponential terms
1342 +        for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1343 +             mol = info_->nextMolecule(mi)) {
1344 +          for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1345 +              atom = mol->nextAtom(ai)) {
1346 +            
1347 +            i = atom->getLocalIndex();
1348 +            if(mm < 0) {
1349 +              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1350 +                     + 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();
1353 +            } else {
1354 +              clm[i] = eCos[l][i].x()*eCos[m][i].y()
1355 +                     - 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();
1358 +            }
1359 +          }
1360 +        }
1361 +        for (int nnn = nMin; nnn <= kLim2; nnn++) {
1362 +          int nn = nnn - kLimit;          
1363 +          int n = abs(nn) + 1;
1364 +          RealType rn = zcl * float(nn);
1365 +          // Test on magnitude of k vector:
1366 +          int kk=ll*ll + mm*mm + nn*nn;
1367 +          if(kk <= kSqLim) {
1368 +            Vector3d kVec = Vector3d(rl, rm, rn);
1369 +            Mat3x3d  k2 = outProduct(kVec, kVec);
1370 +            // Calculate exp(ikr) terms
1371 +            for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1372 +                 mol = info_->nextMolecule(mi)) {
1373 +              for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1374 +                  atom = mol->nextAtom(ai)) {
1375 +                i = atom->getLocalIndex();
1376 +                
1377 +                if (nn < 0) {
1378 +                  ckr[i]=clm[i]*eCos[n][i].z()+slm[i]*eSin[n][i].z();
1379 +                  skr[i]=slm[i]*eCos[n][i].z()-clm[i]*eSin[n][i].z();
1380 +                } else {
1381 +                  ckr[i]=clm[i]*eCos[n][i].z()-slm[i]*eSin[n][i].z();
1382 +                  skr[i]=slm[i]*eCos[n][i].z()+clm[i]*eSin[n][i].z();
1383 +                }
1384 +              }
1385 +            }
1386 +            
1387 +            // Calculate scalar and vector products for each site:
1388 +            
1389 +            for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1390 +                 mol = info_->nextMolecule(mi)) {
1391 +              for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1392 +                  atom = mol->nextAtom(ai)) {
1393 +                i = atom->getGlobalIndex();
1394 +                int atid = atom->getAtomType()->getIdent();
1395 +                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1396 +                              
1397 +                if (data.is_Charge) {
1398 +                  RealType C = data.fixedCharge;
1399 +                  if (atom->isFluctuatingCharge()) C += atom->getFlucQPos();
1400 +                  ckc[i] = C * ckr[i];
1401 +                  cks[i] = C * cks[i];
1402 +                }
1403 +                
1404 +                if (data.is_Dipole) {
1405 +                  Vector3d D = atom->getDipole() * mPoleConverter;
1406 +                  RealType dk = dot(kVec, D);
1407 +                  dxk[i] = cross(kVec, D);
1408 +                  dkc[i] = dk * ckr[i];
1409 +                  dks[i] = dk * skr[i];
1410 +                }
1411 +                if (data.is_Quadrupole) {
1412 +                  Mat3x3d Q = atom->getQuadrupole();
1413 +                  Q *= mPoleConverter;
1414 +                  RealType qk = -( Q * k2 ).trace();
1415 +                  qxk[i] = -2.0 * cross(k2, Q);
1416 +                  qkc[i] = qk * ckr[i];
1417 +                  qks[i] = qk * skr[i];
1418 +                }              
1419 +              }
1420 +            }
1421 +            
1422 +            // calculate vector sums
1423 +            
1424 +            RealType ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1425 +            RealType ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1426 +            RealType dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1427 +            RealType dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1428 +            RealType qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1429 +            RealType qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1430 +            
1431 + #ifdef IS_MPI
1432 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE,
1433 +                                      MPI::SUM);
1434 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckss, 1, MPI::REALTYPE,
1435 +                                      MPI::SUM);
1436 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkcs, 1, MPI::REALTYPE,
1437 +                                      MPI::SUM);
1438 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkss, 1, MPI::REALTYPE,
1439 +                                      MPI::SUM);
1440 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkcs, 1, MPI::REALTYPE,
1441 +                                      MPI::SUM);
1442 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkss, 1, MPI::REALTYPE,
1443 +                                      MPI::SUM);
1444 + #endif        
1445 +            
1446 +            // Accumulate potential energy and virial contribution:
1447 +            
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 +            
1461 +            // Calculate force and torque for each site:
1462 +            
1463 +            for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1464 +                 mol = info_->nextMolecule(mi)) {
1465 +              for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1466 +                  atom = mol->nextAtom(ai)) {
1467 +                
1468 +                i = atom->getLocalIndex();
1469 +                int atid = atom->getAtomType()->getIdent();
1470 +                ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]];
1471 +                
1472 +                RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs)
1473 +                                        - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss));
1474 +                RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
1475 +                                         -ckr[i]*(ckss+dkcs-qkss));
1476 +                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)+
1477 +                                             skr[i]*(ckss+dkcs-qkss));
1478 +                
1479 +                
1480 +                atom->addFrc( 4.0 * rvol * qfrc * kVec );
1481 +                
1482 +                if (data.is_Dipole) {
1483 +                  atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] );
1484 +                }
1485 +                if (data.is_Quadrupole) {
1486 +                  atom->addTrq( 4.0 * rvol * qtrq2 * qxk[i] );
1487 +                }
1488 +              }
1489 +            }
1490 +          }
1491 +        }
1492 +      }
1493 +    }
1494 +  }
1495   }

Diff Legend

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