ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/nonbonded/Electrostatic.cpp
(Generate patch)

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1720 by gezelter, Thu May 24 01:48:29 2012 UTC vs.
Revision 1721 by gezelter, Thu May 24 14:17:42 2012 UTC

# Line 321 | Line 321 | namespace OpenMD {
321        electrostaticAtomData.hardness = fqa.getHardness();
322        electrostaticAtomData.slaterN = fqa.getSlaterN();
323        electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
324 +    } else {
325 +      electrostaticAtomData.is_Fluctuating = false;
326      }
327  
328      pair<map<int,AtomType*>::iterator,bool> ret;    
# Line 364 | Line 366 | namespace OpenMD {
366            rval = RealType(i) * dr;
367            rvals.push_back(rval);
368            J1vals.push_back( sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) );
369 +          // may not be necessary if Slater coulomb integral is symmetric
370            J2vals.push_back( sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) );
371          }
372  
# Line 447 | Line 450 | namespace OpenMD {
450      Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
451  
452      pair<RealType, RealType> res;
453 +    
454 +    // splines for coulomb integrals
455 +    CubicSpline* J1;
456 +    CubicSpline* J2;
457      
458      if (!initialized_) initialize();
459      
# Line 464 | Line 471 | namespace OpenMD {
471      bool i_is_Dipole = data1.is_Dipole;
472      bool i_is_SplitDipole = data1.is_SplitDipole;
473      bool i_is_Quadrupole = data1.is_Quadrupole;
474 +    bool i_is_Fluctuating = data1.is_Fluctuating;
475  
476      bool j_is_Charge = data2.is_Charge;
477      bool j_is_Dipole = data2.is_Dipole;
478      bool j_is_SplitDipole = data2.is_SplitDipole;
479      bool j_is_Quadrupole = data2.is_Quadrupole;
480 +    bool j_is_Fluctuating = data2.is_Fluctuating;
481      
482      if (i_is_Charge) {
483        q_i = data1.fixedCharge;
484 +
485 +      if (i_is_Fluctuating) {
486 +        q_i += *(idat.flucQ1);
487 +      }
488 +      
489        if (idat.excluded) {
490          *(idat.skippedCharge2) += q_i;
491        }
# Line 510 | Line 524 | namespace OpenMD {
524  
525      if (j_is_Charge) {
526        q_j = data2.fixedCharge;
527 +
528 +      if (i_is_Fluctuating)
529 +        q_j += *(idat.flucQ2);
530 +
531        if (idat.excluded) {
532          *(idat.skippedCharge1) += q_j;
533        }
# Line 547 | Line 565 | namespace OpenMD {
565        duduz_j = V3Zero;
566      }
567      
568 +    if (i_is_Fluctuating && j_is_Fluctuating) {
569 +      J1 = Jij[idat.atypes];
570 +      J2 = Jij[make_pair(idat.atypes.second, idat.atypes.first)];
571 +    }
572 +
573      epot = 0.0;
574      dVdr = V3Zero;
575      
# Line 598 | Line 621 | namespace OpenMD {
621  
622            vterm = preVal * riji * erfcVal;          
623            dudr  = -  *(idat.sw)  * preVal * c2;
624 +          
625 +        }
626  
627 +        
628 +        if (i_is_Fluctuating) {
629 +          if (!idat.excluded)
630 +            *(idat.dVdFQ1) += *(idat.sw) * vterm / q_i;
631 +          else {
632 +            res = J1->getValueAndDerivativeAt( *(idat.rij) );
633 +            *(idat.dVdFQ1) += pre11_ * res.first * q_j;
634 +          }
635          }
636 +        if (j_is_Fluctuating) {
637 +          if (!idat.excluded)
638 +            *(idat.dVdFQ2) += *(idat.sw) * vterm / q_j;
639 +          else {
640 +            res = J2->getValueAndDerivativeAt( *(idat.rij) );
641 +            *(idat.dVdFQ2) += pre11_ * res.first * q_i;
642 +          }
643 +        }
644  
645          vpair += vterm;
646          epot +=  *(idat.sw)  * vterm;

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines