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 1601 by gezelter, Thu Aug 4 20:04:35 2011 UTC vs.
Revision 1718 by gezelter, Thu May 24 01:29:59 2012 UTC

# Line 36 | Line 36
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124 234104 (2006).          
38   * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include <stdio.h>
# Line 46 | Line 47
47   #include "nonbonded/Electrostatic.hpp"
48   #include "utils/simError.h"
49   #include "types/NonBondedInteractionType.hpp"
50 < #include "types/DirectionalAtomType.hpp"
50 > #include "types/FixedChargeAdapter.hpp"
51 > #include "types/MultipoleAdapter.hpp"
52   #include "io/Globals.hpp"
53 + #include "nonbonded/SlaterIntegrals.hpp"
54 + #include "utils/PhysicalConstants.hpp"
55  
56 +
57   namespace OpenMD {
58    
59    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
# Line 64 | Line 69 | namespace OpenMD {
69      Globals* simParams_ = info_->getSimParams();
70  
71      summationMap_["HARD"]               = esm_HARD;
72 +    summationMap_["NONE"]               = esm_HARD;
73      summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
74      summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
75      summationMap_["SHIFTED_FORCE"]      = esm_SHIFTED_FORCE;    
# Line 116 | Line 122 | namespace OpenMD {
122          sprintf( painCave.errMsg,
123                   "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
124                   "\t(Input file specified %s .)\n"
125 <                 "\telectrostaticSummationMethod must be one of: \"none\",\n"
125 >                 "\telectrostaticSummationMethod must be one of: \"hard\",\n"
126                   "\t\"shifted_potential\", \"shifted_force\", or \n"
127                   "\t\"reaction_field\".\n", myMethod.c_str() );
128          painCave.isFatal = 1;
# Line 249 | Line 255 | namespace OpenMD {
255        preRF2_ = 2.0 * preRF_;
256      }
257      
258 <    RealType dx = cutoffRadius_ / RealType(np_ - 1);
258 >    // Add a 2 angstrom safety window to deal with cutoffGroups that
259 >    // have charged atoms longer than the cutoffRadius away from each
260 >    // other.  Splining may not be the best choice here.  Direct calls
261 >    // to erfc might be preferrable.
262 >
263 >    RealType dx = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
264      RealType rval;
265      vector<RealType> rvals;
266      vector<RealType> yvals;
# Line 273 | Line 284 | namespace OpenMD {
284      electrostaticAtomData.is_SplitDipole = false;
285      electrostaticAtomData.is_Quadrupole = false;
286  
287 <    if (atomType->isCharge()) {
277 <      GenericData* data = atomType->getPropertyByName("Charge");
287 >    FixedChargeAdapter fca = FixedChargeAdapter(atomType);
288  
289 <      if (data == NULL) {
280 <        sprintf( painCave.errMsg, "Electrostatic::addType could not find "
281 <                 "Charge\n"
282 <                 "\tparameters for atomType %s.\n",
283 <                 atomType->getName().c_str());
284 <        painCave.severity = OPENMD_ERROR;
285 <        painCave.isFatal = 1;
286 <        simError();                  
287 <      }
288 <      
289 <      DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
290 <      if (doubleData == NULL) {
291 <        sprintf( painCave.errMsg,
292 <                 "Electrostatic::addType could not convert GenericData to "
293 <                 "Charge for\n"
294 <                 "\tatom type %s\n", atomType->getName().c_str());
295 <        painCave.severity = OPENMD_ERROR;
296 <        painCave.isFatal = 1;
297 <        simError();          
298 <      }
289 >    if (fca.isFixedCharge()) {
290        electrostaticAtomData.is_Charge = true;
291 <      electrostaticAtomData.charge = doubleData->getData();          
291 >      electrostaticAtomData.charge = fca.getCharge();
292      }
293  
294 <    if (atomType->isDirectional()) {
295 <      DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
296 <      
297 <      if (daType->isDipole()) {
298 <        GenericData* data = daType->getPropertyByName("Dipole");
308 <        
309 <        if (data == NULL) {
310 <          sprintf( painCave.errMsg,
311 <                   "Electrostatic::addType could not find Dipole\n"
312 <                   "\tparameters for atomType %s.\n",
313 <                   daType->getName().c_str());
314 <          painCave.severity = OPENMD_ERROR;
315 <          painCave.isFatal = 1;
316 <          simError();                  
317 <        }
318 <      
319 <        DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
320 <        if (doubleData == NULL) {
321 <          sprintf( painCave.errMsg,
322 <                   "Electrostatic::addType could not convert GenericData to "
323 <                   "Dipole Moment\n"
324 <                   "\tfor atom type %s\n", daType->getName().c_str());
325 <          painCave.severity = OPENMD_ERROR;
326 <          painCave.isFatal = 1;
327 <          simError();          
328 <        }
329 <        electrostaticAtomData.is_Dipole = true;
330 <        electrostaticAtomData.dipole_moment = doubleData->getData();
294 >    MultipoleAdapter ma = MultipoleAdapter(atomType);
295 >    if (ma.isMultipole()) {
296 >      if (ma.isDipole()) {
297 >        electrostaticAtomData.is_Dipole = true;
298 >        electrostaticAtomData.dipole_moment = ma.getDipoleMoment();
299        }
300 <
333 <      if (daType->isSplitDipole()) {
334 <        GenericData* data = daType->getPropertyByName("SplitDipoleDistance");
335 <        
336 <        if (data == NULL) {
337 <          sprintf(painCave.errMsg,
338 <                  "Electrostatic::addType could not find SplitDipoleDistance\n"
339 <                  "\tparameter for atomType %s.\n",
340 <                  daType->getName().c_str());
341 <          painCave.severity = OPENMD_ERROR;
342 <          painCave.isFatal = 1;
343 <          simError();                  
344 <        }
345 <      
346 <        DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
347 <        if (doubleData == NULL) {
348 <          sprintf( painCave.errMsg,
349 <                   "Electrostatic::addType could not convert GenericData to "
350 <                   "SplitDipoleDistance for\n"
351 <                   "\tatom type %s\n", daType->getName().c_str());
352 <          painCave.severity = OPENMD_ERROR;
353 <          painCave.isFatal = 1;
354 <          simError();          
355 <        }
300 >      if (ma.isSplitDipole()) {
301          electrostaticAtomData.is_SplitDipole = true;
302 <        electrostaticAtomData.split_dipole_distance = doubleData->getData();
302 >        electrostaticAtomData.split_dipole_distance = ma.getSplitDipoleDistance();
303        }
304 <
360 <      if (daType->isQuadrupole()) {
361 <        GenericData* data = daType->getPropertyByName("QuadrupoleMoments");
362 <        
363 <        if (data == NULL) {
364 <          sprintf( painCave.errMsg,
365 <                   "Electrostatic::addType could not find QuadrupoleMoments\n"
366 <                   "\tparameter for atomType %s.\n",
367 <                   daType->getName().c_str());
368 <          painCave.severity = OPENMD_ERROR;
369 <          painCave.isFatal = 1;
370 <          simError();                  
371 <        }
372 <        
304 >      if (ma.isQuadrupole()) {
305          // Quadrupoles in OpenMD are set as the diagonal elements
306          // of the diagonalized traceless quadrupole moment tensor.
307          // The column vectors of the unitary matrix that diagonalizes
308          // the quadrupole moment tensor become the eFrame (or the
309          // electrostatic version of the body-fixed frame.
378 –
379 –        Vector3dGenericData* v3dData = dynamic_cast<Vector3dGenericData*>(data);
380 –        if (v3dData == NULL) {
381 –          sprintf( painCave.errMsg,
382 –                   "Electrostatic::addType could not convert GenericData to "
383 –                   "Quadrupole Moments for\n"
384 –                   "\tatom type %s\n", daType->getName().c_str());
385 –          painCave.severity = OPENMD_ERROR;
386 –          painCave.isFatal = 1;
387 –          simError();          
388 –        }
310          electrostaticAtomData.is_Quadrupole = true;
311 <        electrostaticAtomData.quadrupole_moments = v3dData->getData();
311 >        electrostaticAtomData.quadrupole_moments = ma.getQuadrupoleMoments();
312        }
313      }
314      
315 <    AtomTypeProperties atp = atomType->getATP();    
315 >    FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
316 >
317 >    if (fqa.isFluctuatingCharge()) {
318 >      electrostaticAtomData.is_FluctuatingCharge = true;
319 >      electrostaticAtomData.electronegativity = fca.getElectronegativity();
320 >      electrostaticAtomData.hardness = fca.getHardness();
321 >      electrostaticAtomData.slaterN = fca.getSlaterN();
322 >      electrostaticAtomData.slaterZeta = fca.getSlaterZeta();
323 >    }
324  
325      pair<map<int,AtomType*>::iterator,bool> ret;    
326 <    ret = ElectrostaticList.insert( pair<int,AtomType*>(atp.ident, atomType) );
326 >    ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
327 >                                                        atomType) );
328      if (ret.second == false) {
329        sprintf( painCave.errMsg,
330                 "Electrostatic already had a previous entry with ident %d\n",
331 <               atp.ident);
331 >               atomType->getIdent() );
332        painCave.severity = OPENMD_INFO;
333        painCave.isFatal = 0;
334        simError();        
335      }
336      
337 <    ElectrostaticMap[atomType] = electrostaticAtomData;    
337 >    ElectrostaticMap[atomType] = electrostaticAtomData;  
338 >
339 >    // Now, iterate over all known types and add to the mixing map:
340 >    
341 >    map<AtomType*, ElectrostaticAtomData>::iterator it;
342 >    for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
343 >      AtomType* atype2 = (*it).first;
344 >      
345 >      if ((*it).is_FluctuatingCharge && electrostaticAtomData.is_FluctuatingCharge) {
346 >        
347 >        RealType a = electrostaticAtomData.slaterZeta;
348 >        RealType b = (*it).slaterZeta;
349 >        int m = electrostaticAtomData.slaterN;
350 >        int n = (*it).slaterN;
351 >
352 >        // Create the spline of the coulombic integral for s-type
353 >        // Slater orbitals.  Add a 2 angstrom safety window to deal
354 >        // with cutoffGroups that have charged atoms longer than the
355 >        // cutoffRadius away from each other.
356 >
357 >        RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
358 >        vector<RealType> rvals;
359 >        vector<RealType> J1vals;
360 >        vector<RealType> J2vals;
361 >        for (int i = 0; i < np_; i++) {
362 >          rval = RealType(i) * dr;
363 >          rvals.push_back(rval);
364 >          J1vals.push_back( sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) );
365 >          J2vals.push_back( sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) );
366 >        }
367 >
368 >        CubicSpline J1 = new CubicSpline();
369 >        J1->addPoints(rvals, J1vals);
370 >        CubicSpline J2 = new CubicSpline();
371 >        J2->addPoints(rvals, J2vals);
372 >        
373 >        pair<AtomType*, AtomType*> key1, key2;
374 >        key1 = make_pair(atomType, atype2);
375 >        key2 = make_pair(atype2, atomType);
376 >        
377 >        Jij[key1] = J1;
378 >        Jij[key2] = J2;
379 >      }
380 >    }
381 >
382      return;
383    }
384    
# Line 454 | Line 428 | namespace OpenMD {
428      RealType c1, c2, c3, c4;
429      RealType erfcVal(1.0), derfcVal(0.0);
430      RealType BigR;
431 +    RealType two(2.0), three(3.0);
432  
433      Vector3d Q_i, Q_j;
434      Vector3d ux_i, uy_i, uz_i;
# Line 578 | Line 553 | namespace OpenMD {
553        if (j_is_Charge) {
554          if (screeningMethod_ == DAMPED) {
555            // assemble the damping variables
556 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
557 <          erfcVal = res.first;
558 <          derfcVal = res.second;
556 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
557 >          //erfcVal = res.first;
558 >          //derfcVal = res.second;
559 >
560 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
561 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
562 >
563            c1 = erfcVal * riji;
564            c2 = (-derfcVal + c1) * riji;
565          } else {
# Line 610 | Line 589 | namespace OpenMD {
589            if (idat.excluded) {
590              indirect_vpair += preVal * rfVal;
591              indirect_Pot += *(idat.sw) * preVal * rfVal;
592 <            indirect_dVdr += *(idat.sw)  * preVal * 2.0 * rfVal  * riji * rhat;
592 >            indirect_dVdr += *(idat.sw)  * preVal * two * rfVal  * riji * rhat;
593            }
594            
595          } else {
# Line 638 | Line 617 | namespace OpenMD {
617            vpair += vterm;
618            epot +=  *(idat.sw)  * vterm;
619  
620 <          dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
620 >          dVdr +=  -preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
621            duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );  
622  
623            // Even if we excluded this pair from direct interactions,
# Line 667 | Line 646 | namespace OpenMD {
646  
647            if (screeningMethod_ == DAMPED) {
648              // assemble the damping variables
649 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
650 <            erfcVal = res.first;
651 <            derfcVal = res.second;
649 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
650 >            //erfcVal = res.first;
651 >            //derfcVal = res.second;
652 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
653 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
654              c1 = erfcVal * ri;
655              c2 = (-derfcVal + c1) * ri;
656              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 704 | Line 685 | namespace OpenMD {
685            
686          if (screeningMethod_ == DAMPED) {
687            // assemble the damping variables
688 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
689 <          erfcVal = res.first;
690 <          derfcVal = res.second;
688 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
689 >          //erfcVal = res.first;
690 >          //derfcVal = res.second;
691 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
692 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
693            c1 = erfcVal * riji;
694            c2 = (-derfcVal + c1) * riji;
695            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 723 | Line 706 | namespace OpenMD {
706          c2ri = c2 * riji;
707          c3ri = c3 * riji;
708          c4rij = c4 *  *(idat.rij) ;
709 <        rhatdot2 = 2.0 * rhat * c3;
709 >        rhatdot2 = two * rhat * c3;
710          rhatc4 = rhat * c4rij;
711  
712          // calculate the potential
# Line 736 | Line 719 | namespace OpenMD {
719                  
720          // calculate derivatives for the forces and torques
721  
722 <        dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (2.0*cx_j*ux_j + rhat)*c3ri) +
723 <                           qyy_j* (cy2*rhatc4 - (2.0*cy_j*uy_j + rhat)*c3ri) +
724 <                           qzz_j* (cz2*rhatc4 - (2.0*cz_j*uz_j + rhat)*c3ri));
722 >        dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (two*cx_j*ux_j + rhat)*c3ri) +
723 >                           qyy_j* (cy2*rhatc4 - (two*cy_j*uy_j + rhat)*c3ri) +
724 >                           qzz_j* (cz2*rhatc4 - (two*cz_j*uz_j + rhat)*c3ri));
725                            
726          dudux_j += preSw * qxx_j * cx_j * rhatdot2;
727          duduy_j += preSw * qyy_j * cy_j * rhatdot2;
# Line 762 | Line 745 | namespace OpenMD {
745            vpair += vterm;
746            epot +=  *(idat.sw)  * vterm;
747            
748 <          dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
748 >          dVdr += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_ * uz_i);
749            
750            duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
751  
# Line 793 | Line 776 | namespace OpenMD {
776              
777            if (screeningMethod_ == DAMPED) {
778              // assemble the damping variables
779 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
780 <            erfcVal = res.first;
781 <            derfcVal = res.second;
779 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
780 >            //erfcVal = res.first;
781 >            //derfcVal = res.second;
782 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
783 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
784              c1 = erfcVal * ri;
785              c2 = (-derfcVal + c1) * ri;
786              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 838 | Line 823 | namespace OpenMD {
823              
824            a1 = 5.0 * ct_i * ct_j - ct_ij;
825              
826 <          dVdr += preSw * 3.0 * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
826 >          dVdr += preSw * three * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
827  
828 <          duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
829 <          duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
828 >          duduz_i += preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
829 >          duduz_j += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_*uz_i);
830  
831            if (idat.excluded) {
832              indirect_vpair +=  - pref * preRF2_ * ct_ij;
# Line 872 | Line 857 | namespace OpenMD {
857            }
858            if (screeningMethod_ == DAMPED) {
859              // assemble damping variables
860 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
861 <            erfcVal = res.first;
862 <            derfcVal = res.second;
860 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
861 >            //erfcVal = res.first;
862 >            //derfcVal = res.second;
863 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
864 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
865              c1 = erfcVal * ri;
866              c2 = (-derfcVal + c1) * ri;
867              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 923 | Line 910 | namespace OpenMD {
910  
911          if (screeningMethod_ == DAMPED) {
912            // assemble the damping variables
913 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
914 <          erfcVal = res.first;
915 <          derfcVal = res.second;
913 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
914 >          //erfcVal = res.first;
915 >          //derfcVal = res.second;
916 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
917 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
918            c1 = erfcVal * riji;
919            c2 = (-derfcVal + c1) * riji;
920            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 942 | Line 931 | namespace OpenMD {
931          c2ri = c2 * riji;
932          c3ri = c3 * riji;
933          c4rij = c4 *  *(idat.rij) ;
934 <        rhatdot2 = 2.0 * rhat * c3;
934 >        rhatdot2 = two * rhat * c3;
935          rhatc4 = rhat * c4rij;
936  
937          // calculate the potential
# Line 956 | Line 945 | namespace OpenMD {
945  
946          // calculate the derivatives for the forces and torques
947  
948 <        dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (2.0*cx_i*ux_i + rhat)*c3ri) +
949 <                          qyy_i* (cy2*rhatc4 - (2.0*cy_i*uy_i + rhat)*c3ri) +
950 <                          qzz_i* (cz2*rhatc4 - (2.0*cz_i*uz_i + rhat)*c3ri));
948 >        dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (two*cx_i*ux_i + rhat)*c3ri) +
949 >                          qyy_i* (cy2*rhatc4 - (two*cy_i*uy_i + rhat)*c3ri) +
950 >                          qzz_i* (cz2*rhatc4 - (two*cz_i*uz_i + rhat)*c3ri));
951  
952          dudux_i += preSw * qxx_i * cx_i *  rhatdot2;
953          duduy_i += preSw * qyy_i * cy_i *  rhatdot2;
# Line 990 | Line 979 | namespace OpenMD {
979  
980        // only accumulate the forces and torques resulting from the
981        // indirect reaction field terms.
982 +
983        *(idat.vpair) += indirect_vpair;
984        (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
985        *(idat.f1) += indirect_dVdr;

Diff Legend

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