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 1668 by gezelter, Fri Jan 6 19:03:05 2012 UTC vs.
Revision 1720 by gezelter, Thu May 24 01:48:29 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/FluctuatingChargeAdapter.hpp"
52 > #include "types/MultipoleAdapter.hpp"
53   #include "io/Globals.hpp"
54 + #include "nonbonded/SlaterIntegrals.hpp"
55 + #include "utils/PhysicalConstants.hpp"
56  
57 +
58   namespace OpenMD {
59    
60    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
# Line 280 | Line 285 | namespace OpenMD {
285      electrostaticAtomData.is_SplitDipole = false;
286      electrostaticAtomData.is_Quadrupole = false;
287  
288 <    if (atomType->isCharge()) {
284 <      GenericData* data = atomType->getPropertyByName("Charge");
288 >    FixedChargeAdapter fca = FixedChargeAdapter(atomType);
289  
290 <      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 <      }
290 >    if (fca.isFixedCharge()) {
291        electrostaticAtomData.is_Charge = true;
292 <      electrostaticAtomData.charge = doubleData->getData();          
292 >      electrostaticAtomData.fixedCharge = fca.getCharge();
293      }
294  
295 <    if (atomType->isDirectional()) {
296 <      DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
297 <      
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 <        }
295 >    MultipoleAdapter ma = MultipoleAdapter(atomType);
296 >    if (ma.isMultipole()) {
297 >      if (ma.isDipole()) {
298          electrostaticAtomData.is_Dipole = true;
299 <        electrostaticAtomData.dipole_moment = doubleData->getData();
299 >        electrostaticAtomData.dipole_moment = ma.getDipoleMoment();
300        }
301 <
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 <        }
301 >      if (ma.isSplitDipole()) {
302          electrostaticAtomData.is_SplitDipole = true;
303 <        electrostaticAtomData.split_dipole_distance = doubleData->getData();
303 >        electrostaticAtomData.split_dipole_distance = ma.getSplitDipoleDistance();
304        }
305 <
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 <        
305 >      if (ma.isQuadrupole()) {
306          // Quadrupoles in OpenMD are set as the diagonal elements
307          // of the diagonalized traceless quadrupole moment tensor.
308          // The column vectors of the unitary matrix that diagonalizes
309          // the quadrupole moment tensor become the eFrame (or the
310          // 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 –        }
311          electrostaticAtomData.is_Quadrupole = true;
312 <        electrostaticAtomData.quadrupole_moments = v3dData->getData();
312 >        electrostaticAtomData.quadrupole_moments = ma.getQuadrupoleMoments();
313        }
314      }
315      
316 <    AtomTypeProperties atp = atomType->getATP();    
316 >    FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
317  
318 +    if (fqa.isFluctuatingCharge()) {
319 +      electrostaticAtomData.is_Fluctuating = true;
320 +      electrostaticAtomData.electronegativity = fqa.getElectronegativity();
321 +      electrostaticAtomData.hardness = fqa.getHardness();
322 +      electrostaticAtomData.slaterN = fqa.getSlaterN();
323 +      electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
324 +    }
325 +
326      pair<map<int,AtomType*>::iterator,bool> ret;    
327 <    ret = ElectrostaticList.insert( pair<int,AtomType*>(atp.ident, atomType) );
327 >    ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
328 >                                                        atomType) );
329      if (ret.second == false) {
330        sprintf( painCave.errMsg,
331                 "Electrostatic already had a previous entry with ident %d\n",
332 <               atp.ident);
332 >               atomType->getIdent() );
333        painCave.severity = OPENMD_INFO;
334        painCave.isFatal = 0;
335        simError();        
336      }
337      
338 <    ElectrostaticMap[atomType] = electrostaticAtomData;    
338 >    ElectrostaticMap[atomType] = electrostaticAtomData;  
339 >
340 >    // Now, iterate over all known types and add to the mixing map:
341 >    
342 >    map<AtomType*, ElectrostaticAtomData>::iterator it;
343 >    for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
344 >      AtomType* atype2 = (*it).first;
345 >      ElectrostaticAtomData eaData2 = (*it).second;
346 >      if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) {
347 >        
348 >        RealType a = electrostaticAtomData.slaterZeta;
349 >        RealType b = eaData2.slaterZeta;
350 >        int m = electrostaticAtomData.slaterN;
351 >        int n = eaData2.slaterN;
352 >
353 >        // Create the spline of the coulombic integral for s-type
354 >        // Slater orbitals.  Add a 2 angstrom safety window to deal
355 >        // with cutoffGroups that have charged atoms longer than the
356 >        // cutoffRadius away from each other.
357 >
358 >        RealType rval;
359 >        RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
360 >        vector<RealType> rvals;
361 >        vector<RealType> J1vals;
362 >        vector<RealType> J2vals;
363 >        for (int i = 0; i < np_; i++) {
364 >          rval = RealType(i) * dr;
365 >          rvals.push_back(rval);
366 >          J1vals.push_back( sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) );
367 >          J2vals.push_back( sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) );
368 >        }
369 >
370 >        CubicSpline* J1 = new CubicSpline();
371 >        J1->addPoints(rvals, J1vals);
372 >        CubicSpline* J2 = new CubicSpline();
373 >        J2->addPoints(rvals, J2vals);
374 >        
375 >        pair<AtomType*, AtomType*> key1, key2;
376 >        key1 = make_pair(atomType, atype2);
377 >        key2 = make_pair(atype2, atomType);
378 >        
379 >        Jij[key1] = J1;
380 >        Jij[key2] = J2;
381 >      }
382 >    }
383 >
384      return;
385    }
386    
# Line 502 | Line 471 | namespace OpenMD {
471      bool j_is_Quadrupole = data2.is_Quadrupole;
472      
473      if (i_is_Charge) {
474 <      q_i = data1.charge;
474 >      q_i = data1.fixedCharge;
475        if (idat.excluded) {
476          *(idat.skippedCharge2) += q_i;
477        }
# Line 540 | Line 509 | namespace OpenMD {
509      }
510  
511      if (j_is_Charge) {
512 <      q_j = data2.charge;
512 >      q_j = data2.fixedCharge;
513        if (idat.excluded) {
514          *(idat.skippedCharge1) += q_j;
515        }
# Line 1053 | Line 1022 | namespace OpenMD {
1022        }
1023      } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1024        if (i_is_Charge) {        
1025 <        chg1 = data.charge;
1025 >        chg1 = data.fixedCharge;
1026          if (screeningMethod_ == DAMPED) {
1027            self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_;
1028          } else {        

Diff Legend

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