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 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC vs.
Revision 1710 by gezelter, Fri May 18 21:44:02 2012 UTC

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