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 1616 by gezelter, Tue Aug 30 15:45: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 279 | Line 281 | namespace OpenMD {
281      electrostaticAtomData.is_SplitDipole = false;
282      electrostaticAtomData.is_Quadrupole = false;
283  
284 <    if (atomType->isCharge()) {
283 <      GenericData* data = atomType->getPropertyByName("Charge");
284 >    FixedChargeAdapter fca = FixedChargeAdapter(atomType);
285  
286 <      if (data == NULL) {
286 <        sprintf( painCave.errMsg, "Electrostatic::addType could not find "
287 <                 "Charge\n"
288 <                 "\tparameters for atomType %s.\n",
289 <                 atomType->getName().c_str());
290 <        painCave.severity = OPENMD_ERROR;
291 <        painCave.isFatal = 1;
292 <        simError();                  
293 <      }
294 <      
295 <      DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
296 <      if (doubleData == NULL) {
297 <        sprintf( painCave.errMsg,
298 <                 "Electrostatic::addType could not convert GenericData to "
299 <                 "Charge for\n"
300 <                 "\tatom type %s\n", atomType->getName().c_str());
301 <        painCave.severity = OPENMD_ERROR;
302 <        painCave.isFatal = 1;
303 <        simError();          
304 <      }
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 <      
312 <      if (daType->isDipole()) {
313 <        GenericData* data = daType->getPropertyByName("Dipole");
314 <        
315 <        if (data == NULL) {
316 <          sprintf( painCave.errMsg,
317 <                   "Electrostatic::addType could not find Dipole\n"
318 <                   "\tparameters for atomType %s.\n",
319 <                   daType->getName().c_str());
320 <          painCave.severity = OPENMD_ERROR;
321 <          painCave.isFatal = 1;
322 <          simError();                  
323 <        }
324 <      
325 <        DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
326 <        if (doubleData == NULL) {
327 <          sprintf( painCave.errMsg,
328 <                   "Electrostatic::addType could not convert GenericData to "
329 <                   "Dipole Moment\n"
330 <                   "\tfor atom type %s\n", daType->getName().c_str());
331 <          painCave.severity = OPENMD_ERROR;
332 <          painCave.isFatal = 1;
333 <          simError();          
334 <        }
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 <
339 <      if (daType->isSplitDipole()) {
340 <        GenericData* data = daType->getPropertyByName("SplitDipoleDistance");
341 <        
342 <        if (data == NULL) {
343 <          sprintf(painCave.errMsg,
344 <                  "Electrostatic::addType could not find SplitDipoleDistance\n"
345 <                  "\tparameter for atomType %s.\n",
346 <                  daType->getName().c_str());
347 <          painCave.severity = OPENMD_ERROR;
348 <          painCave.isFatal = 1;
349 <          simError();                  
350 <        }
351 <      
352 <        DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
353 <        if (doubleData == NULL) {
354 <          sprintf( painCave.errMsg,
355 <                   "Electrostatic::addType could not convert GenericData to "
356 <                   "SplitDipoleDistance for\n"
357 <                   "\tatom type %s\n", daType->getName().c_str());
358 <          painCave.severity = OPENMD_ERROR;
359 <          painCave.isFatal = 1;
360 <          simError();          
361 <        }
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 <
366 <      if (daType->isQuadrupole()) {
367 <        GenericData* data = daType->getPropertyByName("QuadrupoleMoments");
368 <        
369 <        if (data == NULL) {
370 <          sprintf( painCave.errMsg,
371 <                   "Electrostatic::addType could not find QuadrupoleMoments\n"
372 <                   "\tparameter for atomType %s.\n",
373 <                   daType->getName().c_str());
374 <          painCave.severity = OPENMD_ERROR;
375 <          painCave.isFatal = 1;
376 <          simError();                  
377 <        }
378 <        
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.
384 –
385 –        Vector3dGenericData* v3dData = dynamic_cast<Vector3dGenericData*>(data);
386 –        if (v3dData == NULL) {
387 –          sprintf( painCave.errMsg,
388 –                   "Electrostatic::addType could not convert GenericData to "
389 –                   "Quadrupole Moments for\n"
390 –                   "\tatom type %s\n", daType->getName().c_str());
391 –          painCave.severity = OPENMD_ERROR;
392 –          painCave.isFatal = 1;
393 –          simError();          
394 –        }
307          electrostaticAtomData.is_Quadrupole = true;
308 <        electrostaticAtomData.quadrupole_moments = v3dData->getData();
308 >        electrostaticAtomData.quadrupole_moments = ma.getQuadrupoleMoments();
309        }
310      }
311      
400 –    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 460 | 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 620 | 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 648 | 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 737 | 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 750 | 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 776 | 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 854 | 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 962 | 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 976 | 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;

Diff Legend

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