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

Comparing trunk/src/flucq/FluctuatingChargeConstraints.cpp (file contents):
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 1913 by jmichalk, Wed Jul 24 20:00:51 2013 UTC

# Line 50 | Line 50 | namespace OpenMD {
50   namespace OpenMD {
51  
52    FluctuatingChargeConstraints::FluctuatingChargeConstraints(SimInfo* info) :
53 <    info_(info), hasFlucQ_(false) {
53 >    info_(info), constrainRegions_(false), hasFlucQ_(false), initialized_(false) {
54      
55 <    if (info_->usesFluctuatingCharges()) {
56 <      if (info_->getNFluctuatingCharges() > 0) {
57 <        hasFlucQ_ = true;
55 >  }
56 >
57 >  void FluctuatingChargeConstraints::initialize(){
58 >    if(info_->usesFluctuatingCharges()){
59 >      if(info_->getNFluctuatingCharges() > 0){
60 >        hasFlucQ_ = true;
61        }
62      }
63 +    initialized_ = true;
64    }
65  
66 +  void FluctuatingChargeConstraints::setConstrainRegions(bool cr) {
67 +    constrainRegions_ = cr;
68 +    regionKeys_.clear();
69 +    regionForce_.clear();
70 +    regionCharges_.clear();
71 +
72 +    if (constrainRegions_) {
73 +      SimInfo::MoleculeIterator i;
74 +      Molecule* mol;
75 +      int reg;
76 +      std::set<int> regions;      
77 +
78 +      for (mol = info_->beginMolecule(i); mol != NULL;
79 +           mol = info_->nextMolecule(i)) {
80 +        reg = mol->getRegion();
81 +        if (reg >= 0) regions.insert(reg);
82 +      }
83 +      // resize the keys vector to the largest found value for regions.
84 +      regionKeys_.resize( *(regions.end()) );
85 +      int which = 0;
86 +      for (std::set<int>::iterator r=regions.begin(); r!=regions.end(); ++r) {
87 +        regionKeys_[ (*r) ] = which;
88 +        which++;
89 +      }
90 +      regionForce_.resize( regionKeys_.size() );
91 +      regionCharges_.resize( regionKeys_.size() );
92 +    }
93 +  }
94 +
95 +
96    void FluctuatingChargeConstraints::applyConstraints() {
97 +    if (!initialized_) initialize();
98      if (!hasFlucQ_) return;
64 –    
99      SimInfo::MoleculeIterator i;
100      Molecule::FluctuatingChargeIterator  j;
101      Molecule* mol;
102      Atom* atom;
103      
104 <    RealType totalFrc, totalMolFrc, constrainedFrc;
71 <    
104 >    RealType totalFrc, totalMolFrc, regionFrc, constrainedFrc;
105      // accumulate the total system fluctuating charge forces
106      totalFrc = 0.0;
107 +    if (constrainRegions_) {
108 +      std::fill(regionForce_.begin(), regionForce_.end(), 0.0);
109 +      std::fill(regionCharges_.begin(), regionCharges_.end(), 0);
110 +    }
111  
112      for (mol = info_->beginMolecule(i); mol != NULL;
113           mol = info_->nextMolecule(i)) {
114  
115 +      int region = mol->getRegion();
116 +          
117        for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
118             atom = mol->nextFluctuatingCharge(j)) {
80 –        totalFrc += atom->getFlucQFrc();
81 –      }
119  
120 +        RealType frc = atom->getFlucQFrc();
121 +        totalFrc += frc;
122 +        if (constrainRegions_) {
123 +          if (region >= 0) {        
124 +            regionForce_[regionKeys_[region]] += frc;
125 +            regionCharges_[regionKeys_[region]] += 1;
126 +          }
127 +        }
128 +      }
129      }
130  
131   #ifdef IS_MPI
# Line 87 | Line 133 | namespace OpenMD {
133      // processors:
134      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &totalFrc, 1, MPI::REALTYPE,
135                                MPI::SUM);
136 +
137 +    if (constrainRegions_) {
138 +      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &regionForce_[0],
139 +                                regionForce_.size(), MPI::REALTYPE, MPI::SUM);
140 +      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &regionCharges_[0],
141 +                                regionCharges_.size(), MPI::INT, MPI::SUM);
142 +    }
143 +
144   #endif
145 <
145 >
146      // divide by the total number of fluctuating charges:
147      totalFrc /= info_->getNFluctuatingCharges();
148 +    
149 +    // do the same in the regions:
150 +    if (constrainRegions_) {
151 +      for (int i = 0; i < regionForce_.size(); ++i)  {
152 +        regionForce_[ i ] /= regionCharges_[ i ];
153 +      }
154 +    }
155  
156      for (mol = info_->beginMolecule(i); mol != NULL;
157           mol = info_->nextMolecule(i)) {    
158        
159 <      totalMolFrc = 0.0;
159 >      if (constrainRegions_) {
160 >        int region = mol->getRegion();
161 >        if (region >= 0)
162 >          regionFrc = regionForce_[regionKeys_[region]];
163 >        else
164 >          regionFrc = 0.0;
165 >      } else {
166 >        regionFrc = 0.0;
167 >      }
168  
169 +      totalMolFrc = 0.0;
170 +      
171        // molecular constraints can be done with a second loop.
172        if (mol->constrainTotalCharge()) {
173          for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
# Line 109 | Line 180 | namespace OpenMD {
180        for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
181             atom = mol->nextFluctuatingCharge(j)) {
182          //constrainedFrc = atom->getFlucQFrc() - totalFrc - totalMolFrc;
183 +
184          constrainedFrc = atom->getFlucQFrc() - totalMolFrc;
185 +
186 +        if (constrainRegions_)
187 +          constrainedFrc -= regionFrc;
188 +        
189          atom->setFlucQFrc(constrainedFrc);
190        }      
191      }

Diff Legend

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