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 1613 by gezelter, Thu Aug 18 20:18:19 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 278 | Line 281 | namespace OpenMD {
281      electrostaticAtomData.is_SplitDipole = false;
282      electrostaticAtomData.is_Quadrupole = false;
283  
284 <    if (atomType->isCharge()) {
282 <      GenericData* data = atomType->getPropertyByName("Charge");
284 >    FixedChargeAdapter fca = FixedChargeAdapter(atomType);
285  
286 <      if (data == NULL) {
285 <        sprintf( painCave.errMsg, "Electrostatic::addType could not find "
286 <                 "Charge\n"
287 <                 "\tparameters for atomType %s.\n",
288 <                 atomType->getName().c_str());
289 <        painCave.severity = OPENMD_ERROR;
290 <        painCave.isFatal = 1;
291 <        simError();                  
292 <      }
293 <      
294 <      DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
295 <      if (doubleData == NULL) {
296 <        sprintf( painCave.errMsg,
297 <                 "Electrostatic::addType could not convert GenericData to "
298 <                 "Charge for\n"
299 <                 "\tatom type %s\n", atomType->getName().c_str());
300 <        painCave.severity = OPENMD_ERROR;
301 <        painCave.isFatal = 1;
302 <        simError();          
303 <      }
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 <      
311 <      if (daType->isDipole()) {
312 <        GenericData* data = daType->getPropertyByName("Dipole");
313 <        
314 <        if (data == NULL) {
315 <          sprintf( painCave.errMsg,
316 <                   "Electrostatic::addType could not find Dipole\n"
317 <                   "\tparameters for atomType %s.\n",
318 <                   daType->getName().c_str());
319 <          painCave.severity = OPENMD_ERROR;
320 <          painCave.isFatal = 1;
321 <          simError();                  
322 <        }
323 <      
324 <        DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
325 <        if (doubleData == NULL) {
326 <          sprintf( painCave.errMsg,
327 <                   "Electrostatic::addType could not convert GenericData to "
328 <                   "Dipole Moment\n"
329 <                   "\tfor atom type %s\n", daType->getName().c_str());
330 <          painCave.severity = OPENMD_ERROR;
331 <          painCave.isFatal = 1;
332 <          simError();          
333 <        }
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 <
338 <      if (daType->isSplitDipole()) {
339 <        GenericData* data = daType->getPropertyByName("SplitDipoleDistance");
340 <        
341 <        if (data == NULL) {
342 <          sprintf(painCave.errMsg,
343 <                  "Electrostatic::addType could not find SplitDipoleDistance\n"
344 <                  "\tparameter for atomType %s.\n",
345 <                  daType->getName().c_str());
346 <          painCave.severity = OPENMD_ERROR;
347 <          painCave.isFatal = 1;
348 <          simError();                  
349 <        }
350 <      
351 <        DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
352 <        if (doubleData == NULL) {
353 <          sprintf( painCave.errMsg,
354 <                   "Electrostatic::addType could not convert GenericData to "
355 <                   "SplitDipoleDistance for\n"
356 <                   "\tatom type %s\n", daType->getName().c_str());
357 <          painCave.severity = OPENMD_ERROR;
358 <          painCave.isFatal = 1;
359 <          simError();          
360 <        }
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 <
365 <      if (daType->isQuadrupole()) {
366 <        GenericData* data = daType->getPropertyByName("QuadrupoleMoments");
367 <        
368 <        if (data == NULL) {
369 <          sprintf( painCave.errMsg,
370 <                   "Electrostatic::addType could not find QuadrupoleMoments\n"
371 <                   "\tparameter for atomType %s.\n",
372 <                   daType->getName().c_str());
373 <          painCave.severity = OPENMD_ERROR;
374 <          painCave.isFatal = 1;
375 <          simError();                  
376 <        }
377 <        
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.
383 –
384 –        Vector3dGenericData* v3dData = dynamic_cast<Vector3dGenericData*>(data);
385 –        if (v3dData == NULL) {
386 –          sprintf( painCave.errMsg,
387 –                   "Electrostatic::addType could not convert GenericData to "
388 –                   "Quadrupole Moments for\n"
389 –                   "\tatom type %s\n", daType->getName().c_str());
390 –          painCave.severity = OPENMD_ERROR;
391 –          painCave.isFatal = 1;
392 –          simError();          
393 –        }
307          electrostaticAtomData.is_Quadrupole = true;
308 <        electrostaticAtomData.quadrupole_moments = v3dData->getData();
308 >        electrostaticAtomData.quadrupole_moments = ma.getQuadrupoleMoments();
309        }
310      }
311      
399 –    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 459 | 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 583 | 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 615 | 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 643 | 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 672 | 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 709 | 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 728 | 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 741 | 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 767 | 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 798 | 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 843 | 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 877 | 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 928 | 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 947 | 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 961 | 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 995 | 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