--- branches/development/src/nonbonded/Electrostatic.cpp 2012/05/24 14:17:42 1721 +++ branches/development/src/nonbonded/Electrostatic.cpp 2012/06/07 12:53:46 1750 @@ -193,8 +193,9 @@ namespace OpenMD { // throw warning sprintf( painCave.errMsg, - "Electrostatic::initialize: dampingAlpha was not specified in the input file.\n" - "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", + "Electrostatic::initialize: dampingAlpha was not specified in the\n" + "\tinput file. A default value of %f (1/ang) will be used for the\n" + "\tcutoff of %f (ang).\n", dampingAlpha_, cutoffRadius_); painCave.severity = OPENMD_INFO; painCave.isFatal = 0; @@ -217,7 +218,6 @@ namespace OpenMD { addType(at); } - cutoffRadius2_ = cutoffRadius_ * cutoffRadius_; rcuti_ = 1.0 / cutoffRadius_; rcuti2_ = rcuti_ * rcuti_; @@ -284,6 +284,7 @@ namespace OpenMD { electrostaticAtomData.is_Dipole = false; electrostaticAtomData.is_SplitDipole = false; electrostaticAtomData.is_Quadrupole = false; + electrostaticAtomData.is_Fluctuating = false; FixedChargeAdapter fca = FixedChargeAdapter(atomType); @@ -321,8 +322,6 @@ namespace OpenMD { electrostaticAtomData.hardness = fqa.getHardness(); electrostaticAtomData.slaterN = fqa.getSlaterN(); electrostaticAtomData.slaterZeta = fqa.getSlaterZeta(); - } else { - electrostaticAtomData.is_Fluctuating = false; } pair::iterator,bool> ret; @@ -365,9 +364,9 @@ namespace OpenMD { for (int i = 0; i < np_; i++) { rval = RealType(i) * dr; rvals.push_back(rval); - J1vals.push_back( sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) ); + J1vals.push_back(electrostaticAtomData.hardness * sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) ); // may not be necessary if Slater coulomb integral is symmetric - J2vals.push_back( sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) ); + J2vals.push_back(eaData2.hardness * sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) ); } CubicSpline* J1 = new CubicSpline(); @@ -449,6 +448,7 @@ namespace OpenMD { Vector3d indirect_dVdr(V3Zero); Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero); + RealType coulInt, vFluc1(0.0), vFluc2(0.0); pair res; // splines for coulomb integrals @@ -525,7 +525,7 @@ namespace OpenMD { if (j_is_Charge) { q_j = data2.fixedCharge; - if (i_is_Fluctuating) + if (j_is_Fluctuating) q_j += *(idat.flucQ2); if (idat.excluded) { @@ -592,7 +592,7 @@ namespace OpenMD { c2 = c1 * riji; } - preVal = *(idat.electroMult) * pre11_ * q_i * q_j; + preVal = *(idat.electroMult) * pre11_; if (summationMethod_ == esm_SHIFTED_POTENTIAL) { vterm = preVal * (c1 - c1c_); @@ -623,28 +623,36 @@ namespace OpenMD { dudr = - *(idat.sw) * preVal * c2; } - + vpair += vterm * q_i * q_j; + epot += *(idat.sw) * vterm * q_i * q_j; + dVdr += dudr * rhat * q_i * q_j; + if (i_is_Fluctuating) { - if (!idat.excluded) - *(idat.dVdFQ1) += *(idat.sw) * vterm / q_i; - else { - res = J1->getValueAndDerivativeAt( *(idat.rij) ); - *(idat.dVdFQ1) += pre11_ * res.first * q_j; + if (idat.excluded) { + // vFluc1 is the difference between the direct coulomb integral + // and the normal 1/r-like interaction between point charges. + coulInt = J1->getValueAt( *(idat.rij) ); + vFluc1 = coulInt - (*(idat.sw) * vterm); + } else { + vFluc1 = 0.0; } + *(idat.dVdFQ1) += ( *(idat.sw) * vterm + vFluc1 ) * q_j; } + if (j_is_Fluctuating) { - if (!idat.excluded) - *(idat.dVdFQ2) += *(idat.sw) * vterm / q_j; - else { - res = J2->getValueAndDerivativeAt( *(idat.rij) ); - *(idat.dVdFQ2) += pre11_ * res.first * q_i; + if (idat.excluded) { + // vFluc2 is the difference between the direct coulomb integral + // and the normal 1/r-like interaction between point charges. + coulInt = J2->getValueAt( *(idat.rij) ); + vFluc2 = coulInt - (*(idat.sw) * vterm); + } else { + vFluc2 = 0.0; } + *(idat.dVdFQ2) += ( *(idat.sw) * vterm + vFluc2 ) * q_i; } - - vpair += vterm; - epot += *(idat.sw) * vterm; - dVdr += dudr * rhat; + + } if (j_is_Dipole) { @@ -717,6 +725,9 @@ namespace OpenMD { duduz_j += -preSw * pot_term * rhat; } + if (i_is_Fluctuating) { + *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i; + } } if (j_is_Quadrupole) { @@ -769,6 +780,10 @@ namespace OpenMD { dudux_j += preSw * qxx_j * cx_j * rhatdot2; duduy_j += preSw * qyy_j * cy_j * rhatdot2; duduz_j += preSw * qzz_j * cz_j * rhatdot2; + if (i_is_Fluctuating) { + *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i; + } + } } @@ -844,7 +859,12 @@ namespace OpenMD { // calculate derivatives for the forces and torques dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3); duduz_i += preSw * pot_term * rhat; + } + + if (j_is_Fluctuating) { + *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j; } + } if (j_is_Dipole) { @@ -995,6 +1015,11 @@ namespace OpenMD { dudux_i += preSw * qxx_i * cx_i * rhatdot2; duduy_i += preSw * qyy_i * cy_i * rhatdot2; duduz_i += preSw * qzz_i * cz_i * rhatdot2; + + if (j_is_Fluctuating) { + *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j; + } + } } @@ -1033,13 +1058,11 @@ namespace OpenMD { *(idat.t2) -= cross(uz_j, indirect_duduz_j); } - return; } void Electrostatic::calcSelfCorrection(SelfData &sdat) { - RealType mu1, preVal, chg1, self; - + RealType mu1, preVal, self; if (!initialized_) initialize(); ElectrostaticAtomData data = ElectrostaticMap[sdat.atype]; @@ -1047,6 +1070,14 @@ namespace OpenMD { // logicals bool i_is_Charge = data.is_Charge; bool i_is_Dipole = data.is_Dipole; + bool i_is_Fluctuating = data.is_Fluctuating; + RealType chg1 = data.fixedCharge; + + if (i_is_Fluctuating) { + chg1 += *(sdat.flucQ); + // dVdFQ is really a force, so this is negative the derivative + *(sdat.dVdFQ) -= *(sdat.flucQ) * data.hardness + data.electronegativity; + } if (summationMethod_ == esm_REACTION_FIELD) { if (i_is_Dipole) { @@ -1063,7 +1094,6 @@ namespace OpenMD { } } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) { if (i_is_Charge) { - chg1 = data.fixedCharge; if (screeningMethod_ == DAMPED) { self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_; } else {