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 1710 by gezelter, Fri May 18 21:44:02 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  
54   namespace OpenMD {
# Line 64 | Line 66 | namespace OpenMD {
66      Globals* simParams_ = info_->getSimParams();
67  
68      summationMap_["HARD"]               = esm_HARD;
69 +    summationMap_["NONE"]               = esm_HARD;
70      summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
71      summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
72      summationMap_["SHIFTED_FORCE"]      = esm_SHIFTED_FORCE;    
# Line 116 | Line 119 | namespace OpenMD {
119          sprintf( painCave.errMsg,
120                   "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
121                   "\t(Input file specified %s .)\n"
122 <                 "\telectrostaticSummationMethod must be one of: \"none\",\n"
122 >                 "\telectrostaticSummationMethod must be one of: \"hard\",\n"
123                   "\t\"shifted_potential\", \"shifted_force\", or \n"
124                   "\t\"reaction_field\".\n", myMethod.c_str() );
125          painCave.isFatal = 1;
# Line 249 | Line 252 | namespace OpenMD {
252        preRF2_ = 2.0 * preRF_;
253      }
254      
255 <    RealType dx = cutoffRadius_ / RealType(np_ - 1);
255 >    // Add a 2 angstrom safety window to deal with cutoffGroups that
256 >    // have charged atoms longer than the cutoffRadius away from each
257 >    // other.  Splining may not be the best choice here.  Direct calls
258 >    // to erfc might be preferrable.
259 >
260 >    RealType dx = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
261      RealType rval;
262      vector<RealType> rvals;
263      vector<RealType> yvals;
# Line 273 | Line 281 | namespace OpenMD {
281      electrostaticAtomData.is_SplitDipole = false;
282      electrostaticAtomData.is_Quadrupole = false;
283  
284 <    if (atomType->isCharge()) {
277 <      GenericData* data = atomType->getPropertyByName("Charge");
284 >    FixedChargeAdapter fca = FixedChargeAdapter(atomType);
285  
286 <      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 <      }
286 >    if (fca.isFixedCharge()) {
287        electrostaticAtomData.is_Charge = true;
288 <      electrostaticAtomData.charge = doubleData->getData();          
288 >      electrostaticAtomData.charge = fca.getCharge();
289      }
290  
291 <    if (atomType->isDirectional()) {
292 <      DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
293 <      
306 <      if (daType->isDipole()) {
307 <        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 <        }
291 >    MultipoleAdapter ma = MultipoleAdapter(atomType);
292 >    if (ma.isMultipole()) {
293 >      if (ma.isDipole()) {
294          electrostaticAtomData.is_Dipole = true;
295 <        electrostaticAtomData.dipole_moment = doubleData->getData();
295 >        electrostaticAtomData.dipole_moment = ma.getDipoleMoment();
296        }
297 <
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 <        }
297 >      if (ma.isSplitDipole()) {
298          electrostaticAtomData.is_SplitDipole = true;
299 <        electrostaticAtomData.split_dipole_distance = doubleData->getData();
299 >        electrostaticAtomData.split_dipole_distance = ma.getSplitDipoleDistance();
300        }
301 <
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 <        
301 >      if (ma.isQuadrupole()) {
302          // Quadrupoles in OpenMD are set as the diagonal elements
303          // of the diagonalized traceless quadrupole moment tensor.
304          // The column vectors of the unitary matrix that diagonalizes
305          // the quadrupole moment tensor become the eFrame (or the
306          // 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 –        }
307          electrostaticAtomData.is_Quadrupole = true;
308 <        electrostaticAtomData.quadrupole_moments = v3dData->getData();
308 >        electrostaticAtomData.quadrupole_moments = ma.getQuadrupoleMoments();
309        }
310      }
311      
394 –    AtomTypeProperties atp = atomType->getATP();    
312  
313      pair<map<int,AtomType*>::iterator,bool> ret;    
314 <    ret = ElectrostaticList.insert( pair<int,AtomType*>(atp.ident, atomType) );
314 >    ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
315 >                                                        atomType) );
316      if (ret.second == false) {
317        sprintf( painCave.errMsg,
318                 "Electrostatic already had a previous entry with ident %d\n",
319 <               atp.ident);
319 >               atomType->getIdent() );
320        painCave.severity = OPENMD_INFO;
321        painCave.isFatal = 0;
322        simError();        
# Line 454 | Line 372 | namespace OpenMD {
372      RealType c1, c2, c3, c4;
373      RealType erfcVal(1.0), derfcVal(0.0);
374      RealType BigR;
375 +    RealType two(2.0), three(3.0);
376  
377      Vector3d Q_i, Q_j;
378      Vector3d ux_i, uy_i, uz_i;
# Line 578 | Line 497 | namespace OpenMD {
497        if (j_is_Charge) {
498          if (screeningMethod_ == DAMPED) {
499            // assemble the damping variables
500 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
501 <          erfcVal = res.first;
502 <          derfcVal = res.second;
500 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
501 >          //erfcVal = res.first;
502 >          //derfcVal = res.second;
503 >
504 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
505 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
506 >
507            c1 = erfcVal * riji;
508            c2 = (-derfcVal + c1) * riji;
509          } else {
# Line 610 | Line 533 | namespace OpenMD {
533            if (idat.excluded) {
534              indirect_vpair += preVal * rfVal;
535              indirect_Pot += *(idat.sw) * preVal * rfVal;
536 <            indirect_dVdr += *(idat.sw)  * preVal * 2.0 * rfVal  * riji * rhat;
536 >            indirect_dVdr += *(idat.sw)  * preVal * two * rfVal  * riji * rhat;
537            }
538            
539          } else {
# Line 638 | Line 561 | namespace OpenMD {
561            vpair += vterm;
562            epot +=  *(idat.sw)  * vterm;
563  
564 <          dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
564 >          dVdr +=  -preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
565            duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );  
566  
567            // Even if we excluded this pair from direct interactions,
# Line 667 | Line 590 | namespace OpenMD {
590  
591            if (screeningMethod_ == DAMPED) {
592              // assemble the damping variables
593 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
594 <            erfcVal = res.first;
595 <            derfcVal = res.second;
593 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
594 >            //erfcVal = res.first;
595 >            //derfcVal = res.second;
596 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
597 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
598              c1 = erfcVal * ri;
599              c2 = (-derfcVal + c1) * ri;
600              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 704 | Line 629 | namespace OpenMD {
629            
630          if (screeningMethod_ == DAMPED) {
631            // assemble the damping variables
632 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
633 <          erfcVal = res.first;
634 <          derfcVal = res.second;
632 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
633 >          //erfcVal = res.first;
634 >          //derfcVal = res.second;
635 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
636 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
637            c1 = erfcVal * riji;
638            c2 = (-derfcVal + c1) * riji;
639            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 723 | Line 650 | namespace OpenMD {
650          c2ri = c2 * riji;
651          c3ri = c3 * riji;
652          c4rij = c4 *  *(idat.rij) ;
653 <        rhatdot2 = 2.0 * rhat * c3;
653 >        rhatdot2 = two * rhat * c3;
654          rhatc4 = rhat * c4rij;
655  
656          // calculate the potential
# Line 736 | Line 663 | namespace OpenMD {
663                  
664          // calculate derivatives for the forces and torques
665  
666 <        dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (2.0*cx_j*ux_j + rhat)*c3ri) +
667 <                           qyy_j* (cy2*rhatc4 - (2.0*cy_j*uy_j + rhat)*c3ri) +
668 <                           qzz_j* (cz2*rhatc4 - (2.0*cz_j*uz_j + rhat)*c3ri));
666 >        dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (two*cx_j*ux_j + rhat)*c3ri) +
667 >                           qyy_j* (cy2*rhatc4 - (two*cy_j*uy_j + rhat)*c3ri) +
668 >                           qzz_j* (cz2*rhatc4 - (two*cz_j*uz_j + rhat)*c3ri));
669                            
670          dudux_j += preSw * qxx_j * cx_j * rhatdot2;
671          duduy_j += preSw * qyy_j * cy_j * rhatdot2;
# Line 762 | Line 689 | namespace OpenMD {
689            vpair += vterm;
690            epot +=  *(idat.sw)  * vterm;
691            
692 <          dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
692 >          dVdr += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_ * uz_i);
693            
694            duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
695  
# Line 793 | Line 720 | namespace OpenMD {
720              
721            if (screeningMethod_ == DAMPED) {
722              // assemble the damping variables
723 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
724 <            erfcVal = res.first;
725 <            derfcVal = res.second;
723 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
724 >            //erfcVal = res.first;
725 >            //derfcVal = res.second;
726 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
727 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
728              c1 = erfcVal * ri;
729              c2 = (-derfcVal + c1) * ri;
730              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 838 | Line 767 | namespace OpenMD {
767              
768            a1 = 5.0 * ct_i * ct_j - ct_ij;
769              
770 <          dVdr += preSw * 3.0 * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
770 >          dVdr += preSw * three * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
771  
772 <          duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
773 <          duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
772 >          duduz_i += preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
773 >          duduz_j += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_*uz_i);
774  
775            if (idat.excluded) {
776              indirect_vpair +=  - pref * preRF2_ * ct_ij;
# Line 872 | Line 801 | namespace OpenMD {
801            }
802            if (screeningMethod_ == DAMPED) {
803              // assemble damping variables
804 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
805 <            erfcVal = res.first;
806 <            derfcVal = res.second;
804 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
805 >            //erfcVal = res.first;
806 >            //derfcVal = res.second;
807 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
808 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
809              c1 = erfcVal * ri;
810              c2 = (-derfcVal + c1) * ri;
811              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 923 | Line 854 | namespace OpenMD {
854  
855          if (screeningMethod_ == DAMPED) {
856            // assemble the damping variables
857 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
858 <          erfcVal = res.first;
859 <          derfcVal = res.second;
857 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
858 >          //erfcVal = res.first;
859 >          //derfcVal = res.second;
860 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
861 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
862            c1 = erfcVal * riji;
863            c2 = (-derfcVal + c1) * riji;
864            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 942 | Line 875 | namespace OpenMD {
875          c2ri = c2 * riji;
876          c3ri = c3 * riji;
877          c4rij = c4 *  *(idat.rij) ;
878 <        rhatdot2 = 2.0 * rhat * c3;
878 >        rhatdot2 = two * rhat * c3;
879          rhatc4 = rhat * c4rij;
880  
881          // calculate the potential
# Line 956 | Line 889 | namespace OpenMD {
889  
890          // calculate the derivatives for the forces and torques
891  
892 <        dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (2.0*cx_i*ux_i + rhat)*c3ri) +
893 <                          qyy_i* (cy2*rhatc4 - (2.0*cy_i*uy_i + rhat)*c3ri) +
894 <                          qzz_i* (cz2*rhatc4 - (2.0*cz_i*uz_i + rhat)*c3ri));
892 >        dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (two*cx_i*ux_i + rhat)*c3ri) +
893 >                          qyy_i* (cy2*rhatc4 - (two*cy_i*uy_i + rhat)*c3ri) +
894 >                          qzz_i* (cz2*rhatc4 - (two*cz_i*uz_i + rhat)*c3ri));
895  
896          dudux_i += preSw * qxx_i * cx_i *  rhatdot2;
897          duduy_i += preSw * qyy_i * cy_i *  rhatdot2;
# Line 990 | Line 923 | namespace OpenMD {
923  
924        // only accumulate the forces and torques resulting from the
925        // indirect reaction field terms.
926 +
927        *(idat.vpair) += indirect_vpair;
928        (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
929        *(idat.f1) += indirect_dVdr;

Diff Legend

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