ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/nonbonded/GB.cpp
(Generate patch)

Comparing branches/development/src/nonbonded/GB.cpp (file contents):
Revision 1688 by gezelter, Wed Mar 14 17:56:01 2012 UTC vs.
Revision 1875 by gezelter, Fri May 17 14:41:42 2013 UTC

# Line 35 | Line 35
35   *                                                                      
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).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39   * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
# Line 46 | Line 46
46   #include <cmath>
47   #include "nonbonded/GB.hpp"
48   #include "utils/simError.h"
49 + #include "types/LennardJonesAdapter.hpp"
50 + #include "types/GayBerneAdapter.hpp"
51  
52   using namespace std;
53   namespace OpenMD {
# Line 87 | Line 89 | namespace OpenMD {
89  
90  
91    GB::GB() : name_("GB"), initialized_(false), mu_(2.0), nu_(1.0), forceField_(NULL) {}
90 –  
91 –  GayBerneParam GB::getGayBerneParam(AtomType* atomType) {
92      
93 –    // Do sanity checking on the AtomType we were passed before
94 –    // building any data structures:
95 –    if (!atomType->isGayBerne()) {
96 –      sprintf( painCave.errMsg,
97 –               "GB::getGayBerneParam was passed an atomType (%s) that does\n"
98 –               "\tnot appear to be a Gay-Berne atom.\n",
99 –               atomType->getName().c_str());
100 –      painCave.severity = OPENMD_ERROR;
101 –      painCave.isFatal = 1;
102 –      simError();
103 –    }
104 –    
105 –    DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
106 –    GenericData* data = daType->getPropertyByName("GayBerne");
107 –    if (data == NULL) {
108 –      sprintf( painCave.errMsg, "GB::getGayBerneParam could not find\n"
109 –               "\tGay-Berne parameters for atomType %s.\n",
110 –               daType->getName().c_str());
111 –      painCave.severity = OPENMD_ERROR;
112 –      painCave.isFatal = 1;
113 –      simError();
114 –    }
115 –    
116 –    GayBerneParamGenericData* gbData = dynamic_cast<GayBerneParamGenericData*>(data);
117 –    if (gbData == NULL) {
118 –      sprintf( painCave.errMsg,
119 –               "GB::getGayBerneParam could not convert GenericData to\n"
120 –               "\tGayBerneParamGenericData for atom type %s\n",
121 –               daType->getName().c_str());
122 –      painCave.severity = OPENMD_ERROR;
123 –      painCave.isFatal = 1;
124 –      simError();          
125 –    }
126 –    
127 –    return gbData->getData();
128 –  }
129 –
130 –  LJParam GB::getLJParam(AtomType* atomType) {
131 –    
132 –    // Do sanity checking on the AtomType we were passed before
133 –    // building any data structures:
134 –    if (!atomType->isLennardJones()) {
135 –      sprintf( painCave.errMsg,
136 –               "GB::getLJParam was passed an atomType (%s) that does not\n"
137 –               "\tappear to be a Lennard-Jones atom.\n",
138 –               atomType->getName().c_str());
139 –      painCave.severity = OPENMD_ERROR;
140 –      painCave.isFatal = 1;
141 –      simError();
142 –    }
143 –    
144 –    GenericData* data = atomType->getPropertyByName("LennardJones");
145 –    if (data == NULL) {
146 –      sprintf( painCave.errMsg, "GB::getLJParam could not find Lennard-Jones\n"
147 –               "\tparameters for atomType %s.\n", atomType->getName().c_str());
148 –      painCave.severity = OPENMD_ERROR;
149 –      painCave.isFatal = 1;
150 –      simError();
151 –    }
152 –    
153 –    LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
154 –    if (ljData == NULL) {
155 –      sprintf( painCave.errMsg,
156 –               "GB::getLJParam could not convert GenericData to LJParam for\n"
157 –               "\tatom type %s\n", atomType->getName().c_str());
158 –      painCave.severity = OPENMD_ERROR;
159 –      painCave.isFatal = 1;
160 –      simError();          
161 –    }
162 –    
163 –    return ljData->getData();
164 –  }
165 –  
166 –  RealType GB::getLJEpsilon(AtomType* atomType) {    
167 –    LJParam ljParam = getLJParam(atomType);
168 –    return ljParam.epsilon;
169 –  }
170 –  RealType GB::getLJSigma(AtomType* atomType) {    
171 –    LJParam ljParam = getLJParam(atomType);
172 –    return ljParam.sigma;
173 –  }
174 –  
93    void GB::initialize() {    
94      
95      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
# Line 187 | Line 105 | namespace OpenMD {
105      for (at = atomTypes->beginType(i); at != NULL;
106           at = atomTypes->nextType(i)) {
107        
108 <      if (at->isGayBerne() || at->isLennardJones())
108 >      LennardJonesAdapter lja = LennardJonesAdapter(at);
109 >      GayBerneAdapter gba = GayBerneAdapter(at);
110 >
111 >      if (gba.isGayBerne() || lja.isLennardJones())
112          addType(at);
113      }
114    
# Line 196 | Line 117 | namespace OpenMD {
117        
118    void GB::addType(AtomType* atomType){
119      // add it to the map:
199 –    AtomTypeProperties atp = atomType->getATP();    
120  
121      pair<map<int,AtomType*>::iterator,bool> ret;    
122 <    ret = GBMap.insert( pair<int, AtomType*>(atp.ident, atomType) );
122 >    ret = GBMap.insert( pair<int, AtomType*>(atomType->getIdent(), atomType) );
123      if (ret.second == false) {
124        sprintf( painCave.errMsg,
125                 "GB already had a previous entry with ident %d\n",
126 <               atp.ident);
126 >               atomType->getIdent() );
127        painCave.severity = OPENMD_INFO;
128        painCave.isFatal = 0;
129        simError();        
130      }
131      
132      RealType d1, l1, eX1, eS1, eE1, dw1;
133 <    
134 <    if (atomType->isGayBerne()) {
135 <      GayBerneParam gb1 = getGayBerneParam(atomType);
136 <      d1 = gb1.GB_d;
137 <      l1 = gb1.GB_l;
138 <      eX1 = gb1.GB_eps_X;
139 <      eS1 = gb1.GB_eps_S;
140 <      eE1 = gb1.GB_eps_E;
141 <      dw1 = gb1.GB_dw;
142 <    } else if (atomType->isLennardJones()) {
143 <      d1 = getLJSigma(atomType) / sqrt(2.0);
133 >        
134 >    LennardJonesAdapter lja1 = LennardJonesAdapter(atomType);
135 >    GayBerneAdapter gba1 = GayBerneAdapter(atomType);
136 >    if (gba1.isGayBerne()) {
137 >      d1 = gba1.getD();
138 >      l1 = gba1.getL();
139 >      eX1 = gba1.getEpsX();
140 >      eS1 = gba1.getEpsS();
141 >      eE1 = gba1.getEpsE();
142 >      dw1 = gba1.getDw();
143 >    } else if (lja1.isLennardJones()) {
144 >      d1 = lja1.getSigma() / sqrt(2.0);
145        l1 = d1;
146 <      eX1 = getLJEpsilon(atomType);
146 >      eX1 = lja1.getEpsilon();
147        eS1 = eX1;
148        eE1 = eX1;
149        dw1 = 1.0;      
# Line 243 | Line 164 | namespace OpenMD {
164      for( it = GBMap.begin(); it != GBMap.end(); ++it) {
165        
166        AtomType* atype2 = (*it).second;
167 <      
167 >      LennardJonesAdapter lja2 = LennardJonesAdapter(atype2);
168 >      GayBerneAdapter gba2 = GayBerneAdapter(atype2);
169        RealType d2, l2, eX2, eS2, eE2, dw2;
170        
171 <      if (atype2->isGayBerne()) {
172 <        GayBerneParam gb2 = getGayBerneParam(atype2);
173 <        d2 = gb2.GB_d;
174 <        l2 = gb2.GB_l;
175 <        eX2 = gb2.GB_eps_X;
176 <        eS2 = gb2.GB_eps_S;
177 <        eE2 = gb2.GB_eps_E;
178 <        dw2 = gb2.GB_dw;
179 <      } else if (atype2->isLennardJones()) {
258 <        d2 = getLJSigma(atype2) / sqrt(2.0);
171 >      if (gba2.isGayBerne()) {
172 >        d2 = gba2.getD();
173 >        l2 = gba2.getL();
174 >        eX2 = gba2.getEpsX();
175 >        eS2 = gba2.getEpsS();
176 >        eE2 = gba2.getEpsE();
177 >        dw2 = gba2.getDw();
178 >      } else if (lja2.isLennardJones()) {
179 >        d2 = lja2.getSigma() / sqrt(2.0);
180          l2 = d2;
181 <        eX2 = getLJEpsilon(atype2);
181 >        eX2 = lja2.getEpsilon();
182          eS2 = eX2;
183          eE2 = eX2;
184          dw2 = 1.0;
185 <      }
185 >      } else {
186 >        sprintf( painCave.errMsg,
187 >                 "GB::addType found an atomType (%s) that does not\n"
188 >                 "\tappear to be a Gay-Berne or Lennard-Jones atom.\n",
189 >                 atype2->getName().c_str());
190 >        painCave.severity = OPENMD_ERROR;
191 >        painCave.isFatal = 1;
192 >        simError();
193 >      }
194 >
195                        
196        GBInteractionData mixer1, mixer2;    
197        
# Line 304 | Line 234 | namespace OpenMD {
234  
235        // only add this pairing if at least one of the atoms is a Gay-Berne atom
236  
237 <      if (atomType->isGayBerne() || atype2->isGayBerne()) {
237 >      if (gba1.isGayBerne() || gba2.isGayBerne()) {
238  
239          pair<AtomType*, AtomType*> key1, key2;
240          key1 = make_pair(atomType, atype2);
# Line 353 | Line 283 | namespace OpenMD {
283  
284      RealType a, b, g;
285  
286 +    // This is not good.  We should store this in the mixing map, and not
287 +    // query atom types in calc force.
288      bool i_is_LJ = idat.atypes.first->isLennardJones();
289      bool j_is_LJ = idat.atypes.second->isLennardJones();
290      
# Line 469 | Line 401 | namespace OpenMD {
401  
402      RealType cut = 0.0;
403  
404 <    if (atypes.first->isGayBerne()) {
405 <      GayBerneParam gb1 = getGayBerneParam(atypes.first);
406 <      RealType d1 = gb1.GB_d;
407 <      RealType l1 = gb1.GB_l;
404 >    LennardJonesAdapter lja1 = LennardJonesAdapter(atypes.first);
405 >    GayBerneAdapter gba1 = GayBerneAdapter(atypes.first);
406 >    LennardJonesAdapter lja2 = LennardJonesAdapter(atypes.second);
407 >    GayBerneAdapter gba2 = GayBerneAdapter(atypes.second);
408 >
409 >    if (gba1.isGayBerne()) {
410 >      RealType d1 = gba1.getD();
411 >      RealType l1 = gba1.getL();
412        // sigma is actually sqrt(2)*l  for prolate ellipsoids
413        cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d1, l1));
414 <    } else if (atypes.first->isLennardJones()) {
415 <      cut = max(cut, RealType(2.5) * getLJSigma(atypes.first));
414 >    } else if (lja1.isLennardJones()) {
415 >      cut = max(cut, RealType(2.5) * lja1.getSigma());
416      }
417  
418 <    if (atypes.second->isGayBerne()) {
419 <      GayBerneParam gb2 = getGayBerneParam(atypes.second);
420 <      RealType d2 = gb2.GB_d;
485 <      RealType l2 = gb2.GB_l;
418 >    if (gba2.isGayBerne()) {
419 >      RealType d2 = gba2.getD();
420 >      RealType l2 = gba2.getL();
421        cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d2, l2));
422 <    } else if (atypes.second->isLennardJones()) {
423 <      cut = max(cut, RealType(2.5) * getLJSigma(atypes.second));
422 >    } else if (lja2.isLennardJones()) {
423 >      cut = max(cut, RealType(2.5) * lja2.getSigma());
424      }
425    
426      return cut;

Diff Legend

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