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 1722 by gezelter, Thu May 24 14:17:42 2012 UTC vs.
Revision 1723 by gezelter, Thu May 24 20:59:54 2012 UTC

# Line 217 | Line 217 | namespace OpenMD {
217          addType(at);
218      }
219      
220 –
220      cutoffRadius2_ = cutoffRadius_ * cutoffRadius_;
221      rcuti_ = 1.0 / cutoffRadius_;
222      rcuti2_ = rcuti_ * rcuti_;
# Line 284 | Line 283 | namespace OpenMD {
283      electrostaticAtomData.is_Dipole = false;
284      electrostaticAtomData.is_SplitDipole = false;
285      electrostaticAtomData.is_Quadrupole = false;
286 +    electrostaticAtomData.is_Fluctuating = false;
287  
288      FixedChargeAdapter fca = FixedChargeAdapter(atomType);
289  
# 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;
324      }
325  
326      pair<map<int,AtomType*>::iterator,bool> ret;    
# Line 449 | Line 447 | namespace OpenMD {
447      Vector3d indirect_dVdr(V3Zero);
448      Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
449  
450 +    RealType coulInt, vFluc1(0.0), vFluc2(0.0);
451      pair<RealType, RealType> res;
452      
453      // splines for coulomb integrals
# Line 623 | Line 622 | namespace OpenMD {
622            dudr  = -  *(idat.sw)  * preVal * c2;
623            
624          }
626 –
625          
626 +        vpair += vterm;
627 +        epot +=  *(idat.sw)  * vterm;
628 +        dVdr += dudr * rhat;
629 +
630          if (i_is_Fluctuating) {
631 <          if (!idat.excluded)
632 <            *(idat.dVdFQ1) += *(idat.sw) * vterm / q_i;
633 <          else {
634 <            res = J1->getValueAndDerivativeAt( *(idat.rij) );
635 <            *(idat.dVdFQ1) += pre11_ * res.first * q_j;
631 >          if (idat.excluded) {
632 >            // vFluc1 is the difference between the direct coulomb integral
633 >            // and the normal 1/r-like  interaction between point charges.
634 >            coulInt = J1->getValueAt( *(idat.rij) );
635 >            vFluc1 = pre11_ * coulInt * q_i * q_j  - (*(idat.sw) * vterm);
636 >          } else {
637 >            vFluc1 = 0.0;
638            }
639 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm + vFluc1 ) / q_i;
640          }
641 +
642          if (j_is_Fluctuating) {
643 <          if (!idat.excluded)
644 <            *(idat.dVdFQ2) += *(idat.sw) * vterm / q_j;
645 <          else {
646 <            res = J2->getValueAndDerivativeAt( *(idat.rij) );
647 <            *(idat.dVdFQ2) += pre11_ * res.first * q_i;
643 >          if (idat.excluded) {
644 >            // vFluc2 is the difference between the direct coulomb integral
645 >            // and the normal 1/r-like  interaction between point charges.
646 >            coulInt = J2->getValueAt( *(idat.rij) );
647 >            vFluc2 = pre11_ * coulInt * q_i * q_j  - (*(idat.sw) * vterm);
648 >          } else {
649 >            vFluc2 = 0.0;
650            }
651 +          *(idat.dVdFQ2) += ( *(idat.sw) * vterm + vFluc2 ) / q_j;
652          }
653 <
654 <        vpair += vterm;
646 <        epot +=  *(idat.sw)  * vterm;
647 <        dVdr += dudr * rhat;                
653 >          
654 >
655        }
656  
657        if (j_is_Dipole) {
# Line 716 | Line 723 | namespace OpenMD {
723            dVdr += -preSw * (uz_j * c2ri - ct_j * rhat * sc2 * c3);
724            duduz_j += -preSw * pot_term * rhat;
725  
726 +        }
727 +        if (i_is_Fluctuating) {
728 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
729          }
730        }
731  
# Line 769 | Line 779 | namespace OpenMD {
779          dudux_j += preSw * qxx_j * cx_j * rhatdot2;
780          duduy_j += preSw * qyy_j * cy_j * rhatdot2;
781          duduz_j += preSw * qzz_j * cz_j * rhatdot2;
782 +        if (i_is_Fluctuating) {
783 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
784 +        }
785 +
786        }
787      }
788      
# Line 845 | Line 859 | namespace OpenMD {
859            dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
860            duduz_i += preSw * pot_term * rhat;
861          }
862 +
863 +        if (j_is_Fluctuating) {
864 +          *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
865 +        }
866 +
867        }
868  
869        if (j_is_Dipole) {
# Line 995 | Line 1014 | namespace OpenMD {
1014          dudux_i += preSw * qxx_i * cx_i *  rhatdot2;
1015          duduy_i += preSw * qyy_i * cy_i *  rhatdot2;
1016          duduz_i += preSw * qzz_i * cz_i *  rhatdot2;
1017 +
1018 +        if (j_is_Fluctuating) {
1019 +          *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
1020 +        }
1021 +
1022        }
1023      }
1024  
# Line 1033 | Line 1057 | namespace OpenMD {
1057          *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1058      }
1059  
1036 –
1060      return;
1061    }  
1062      

Diff Legend

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