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

Comparing trunk/src/nonbonded/RepulsivePower.cpp (file contents):
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 1895 by gezelter, Mon Jul 1 21:09:37 2013 UTC

# Line 57 | Line 57 | namespace OpenMD {
57    
58    void RepulsivePower::initialize() {    
59  
60 +    RPtypes.clear();
61 +    RPtids.clear();
62 +    MixingMap.clear();
63 +    RPtids.resize( forceField_->getNAtomType(), -1);
64 +
65      ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
66      ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
67      ForceField::NonBondedInteractionTypeContainer::KeyType keys;
68      NonBondedInteractionType* nbt;
69 +    int rptid1, rptid2;
70  
71      for (nbt = nbiTypes->beginType(j); nbt != NULL;
72           nbt = nbiTypes->nextType(j)) {
# Line 70 | Line 76 | namespace OpenMD {
76          AtomType* at1 = forceField_->getAtomType(keys[0]);
77          AtomType* at2 = forceField_->getAtomType(keys[1]);
78  
79 <
79 >        int atid1 = at1->getIdent();
80 >        if (RPtids[atid1] == -1) {
81 >          rptid1 = RPtypes.size();
82 >          RPtypes.insert(atid1);
83 >          RPtids[atid1] = rptid1;        
84 >        }
85 >        int atid2 = at2->getIdent();
86 >        if (RPtids[atid2] == -1) {
87 >          rptid2 = RPtypes.size();
88 >          RPtypes.insert(atid2);
89 >          RPtids[atid2] = rptid2;        
90 >        }
91 >          
92          RepulsivePowerInteractionType* rpit = dynamic_cast<RepulsivePowerInteractionType*>(nbt);
75 –
93          if (rpit == NULL) {
94            sprintf( painCave.errMsg,
95                     "RepulsivePower::initialize could not convert NonBondedInteractionType\n"
# Line 83 | Line 100 | namespace OpenMD {
100            painCave.isFatal = 1;
101            simError();          
102          }
86 –
103          
104          RealType sigma = rpit->getSigma();
105          RealType epsilon = rpit->getEpsilon();
# Line 107 | Line 123 | namespace OpenMD {
123      mixer.sigmai = 1.0 / mixer.sigma;
124      mixer.nRep = nRep;
125  
126 <    pair<AtomType*, AtomType*> key1, key2;
127 <    key1 = make_pair(atype1, atype2);
128 <    key2 = make_pair(atype2, atype1);
126 >    int rptid1 = RPtids[atype1->getIdent()];
127 >    int rptid2 = RPtids[atype2->getIdent()];
128 >    int nRP = RPtypes.size();
129 >
130 >    MixingMap.resize(nRP);
131 >    MixingMap[rptid1].resize(nRP);
132      
133 <    MixingMap[key1] = mixer;
134 <    if (key2 != key1) {
135 <      MixingMap[key2] = mixer;
133 >    MixingMap[rptid1][rptid2] = mixer;
134 >    if (rptid2 != rptid1) {
135 >      MixingMap[rptid2].resize(nRP);
136 >      MixingMap[rptid2][rptid1] = mixer;
137      }    
138    }
139    
140    void RepulsivePower::calcForce(InteractionData &idat) {
141  
142      if (!initialized_) initialize();
123 –    
124 –    map<pair<AtomType*, AtomType*>, RPInteractionData>::iterator it;
125 –    it = MixingMap.find( idat.atypes );
143  
144 <    if (it != MixingMap.end()) {
145 <
146 <      RPInteractionData mixer = (*it).second;
147 <      RealType sigmai = mixer.sigmai;
148 <      RealType epsilon = mixer.epsilon;
149 <      int nRep = mixer.nRep;
150 <      
151 <      RealType ros;
152 <      RealType rcos;
153 <      RealType myPot = 0.0;
154 <      RealType myPotC = 0.0;
155 <      RealType myDeriv = 0.0;
156 <      RealType myDerivC = 0.0;
157 <    
158 <      ros = *(idat.rij) * sigmai;    
159 <      
160 <      getNRepulsionFunc(ros, nRep, myPot, myDeriv);
161 <      
162 <      if (idat.shiftedPot) {
163 <        rcos = *(idat.rcut) * sigmai;
164 <        getNRepulsionFunc(rcos, nRep, myPotC, myDerivC);
165 <        myDerivC = 0.0;
166 <      } else if (idat.shiftedForce) {
167 <        rcos = *(idat.rcut) * sigmai;
168 <        getNRepulsionFunc(rcos, nRep, myPotC, myDerivC);
169 <        myPotC = myPotC + myDerivC * (*(idat.rij) - *(idat.rcut)) * sigmai;
170 <      } else {
154 <        myPotC = 0.0;
155 <        myDerivC = 0.0;        
156 <      }
157 <
158 <      RealType pot_temp = *(idat.vdwMult) * epsilon * (myPot - myPotC);
159 <      *(idat.vpair) += pot_temp;
160 <      
161 <      RealType dudr = *(idat.sw) * *(idat.vdwMult) * epsilon * (myDeriv -
162 <                                                                myDerivC)*sigmai;      
163 <
164 <      (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
165 <      *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
144 >    RPInteractionData &mixer = MixingMap[RPtids[idat.atid1]][RPtids[idat.atid2]];
145 >    RealType sigmai = mixer.sigmai;
146 >    RealType epsilon = mixer.epsilon;
147 >    int nRep = mixer.nRep;
148 >    
149 >    RealType ros;
150 >    RealType rcos;
151 >    RealType myPot = 0.0;
152 >    RealType myPotC = 0.0;
153 >    RealType myDeriv = 0.0;
154 >    RealType myDerivC = 0.0;
155 >    
156 >    ros = *(idat.rij) * sigmai;    
157 >    
158 >    getNRepulsionFunc(ros, nRep, myPot, myDeriv);
159 >    
160 >    if (idat.shiftedPot) {
161 >      rcos = *(idat.rcut) * sigmai;
162 >      getNRepulsionFunc(rcos, nRep, myPotC, myDerivC);
163 >      myDerivC = 0.0;
164 >    } else if (idat.shiftedForce) {
165 >      rcos = *(idat.rcut) * sigmai;
166 >      getNRepulsionFunc(rcos, nRep, myPotC, myDerivC);
167 >      myPotC = myPotC + myDerivC * (*(idat.rij) - *(idat.rcut)) * sigmai;
168 >    } else {
169 >      myPotC = 0.0;
170 >      myDerivC = 0.0;        
171      }
172 +    
173 +    RealType pot_temp = *(idat.vdwMult) * epsilon * (myPot - myPotC);
174 +    *(idat.vpair) += pot_temp;
175 +    
176 +    RealType dudr = *(idat.sw) * *(idat.vdwMult) * epsilon * (myDeriv -
177 +                                                              myDerivC)*sigmai;      
178 +    
179 +    (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
180 +    *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
181 +    
182      return;
183    }
184  
185 <  void RepulsivePower::getNRepulsionFunc(RealType r, int n, RealType &pot, RealType &deriv) {
185 >  void RepulsivePower::getNRepulsionFunc(const RealType &r, int &n, RealType &pot, RealType &deriv) {
186  
187      RealType ri = 1.0 / r;
188      RealType rin = pow(ri, n);
# Line 182 | Line 197 | namespace OpenMD {
197  
198    RealType RepulsivePower::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {  
199      if (!initialized_) initialize();  
200 <    map<pair<AtomType*, AtomType*>, RPInteractionData>::iterator it;
201 <    it = MixingMap.find(atypes);
202 <    if (it == MixingMap.end())
203 <      return 0.0;
204 <    else  {
205 <      RPInteractionData mixer = (*it).second;
200 >    
201 >    int atid1 = atypes.first->getIdent();
202 >    int atid2 = atypes.second->getIdent();
203 >    int rptid1 = RPtids[atid1];
204 >    int rptid2 = RPtids[atid2];
205 >    
206 >    if (rptid1 == -1 || rptid2 == -1) return 0.0;
207 >    else {      
208 >      RPInteractionData mixer = MixingMap[rptid1][rptid2];
209        return 2.5 * mixer.sigma;
210      }
211 <  }
194 <
211 >  }  
212   }
213  

Diff Legend

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