--- branches/development/src/nonbonded/Electrostatic.cpp 2012/05/24 14:17:42 1721 +++ branches/development/src/nonbonded/Electrostatic.cpp 2012/06/22 20:01:37 1761 @@ -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; @@ -362,12 +361,13 @@ namespace OpenMD { vector rvals; vector J1vals; vector J2vals; - for (int i = 0; i < np_; i++) { + // don't start at i = 0, as rval = 0 is undefined for the slater overlap integrals. + for (int i = 1; i < np_+1; 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 +449,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 +526,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) { @@ -578,12 +579,12 @@ namespace OpenMD { if (j_is_Charge) { if (screeningMethod_ == DAMPED) { // assemble the damping variables - //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) ); - //erfcVal = res.first; - //derfcVal = res.second; + res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) ); + erfcVal = res.first; + derfcVal = res.second; - erfcVal = erfc(dampingAlpha_ * *(idat.rij)); - derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2)); + //erfcVal = erfc(dampingAlpha_ * *(idat.rij)); + //derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2)); c1 = erfcVal * riji; c2 = (-derfcVal + c1) * riji; @@ -592,7 +593,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 +624,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 +726,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 +781,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; + } + } } @@ -845,6 +861,11 @@ namespace OpenMD { 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 +1016,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; + } + } } @@ -1024,6 +1050,9 @@ namespace OpenMD { // indirect reaction field terms. *(idat.vpair) += indirect_vpair; + + (*(idat.excludedPot))[ELECTROSTATIC_FAMILY] += (*(idat.sw) * vterm + + vFluc1 ) * q_i * q_j; (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot; *(idat.f1) += indirect_dVdr; @@ -1033,13 +1062,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 +1074,16 @@ 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; + (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) * + (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity); + } if (summationMethod_ == esm_REACTION_FIELD) { if (i_is_Dipole) { @@ -1063,7 +1100,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 {