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

Comparing trunk/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1895 by gezelter, Mon Jul 1 21:09:37 2013 UTC vs.
Revision 1932 by gezelter, Wed Aug 21 16:52:56 2013 UTC

# Line 44 | Line 44
44   #include <string.h>
45  
46   #include <cmath>
47 + #include <numeric>
48   #include "nonbonded/Electrostatic.hpp"
49   #include "utils/simError.h"
50   #include "types/NonBondedInteractionType.hpp"
# Line 55 | Line 56
56   #include "utils/PhysicalConstants.hpp"
57   #include "math/erfc.hpp"
58   #include "math/SquareMatrix.hpp"
59 + #include "primitives/Molecule.hpp"
60 + #ifdef IS_MPI
61 + #include <mpi.h>
62 + #endif
63  
64   namespace OpenMD {
65    
# Line 191 | Line 196 | namespace OpenMD {
196        simError();
197      }
198            
199 <    if (screeningMethod_ == DAMPED) {      
199 >    if (screeningMethod_ == DAMPED || summationMethod_ == esm_EWALD_FULL) {
200        if (!simParams_->haveDampingAlpha()) {
201          // first set a cutoff dependent alpha value
202          // we assume alpha depends linearly with rcut from 0 to 20.5 ang
203          dampingAlpha_ = 0.425 - cutoffRadius_* 0.02;
204 <        if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;
200 <        
204 >        if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;        
205          // throw warning
206          sprintf( painCave.errMsg,
207                   "Electrostatic::initialize: dampingAlpha was not specified in the\n"
# Line 213 | Line 217 | namespace OpenMD {
217        haveDampingAlpha_ = true;
218      }
219  
220 +
221      Etypes.clear();
222      Etids.clear();
223      FQtypes.clear();
# Line 262 | Line 267 | namespace OpenMD {
267        b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2;
268        b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2;
269        b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2;
270 <      selfMult_ = b0c + a2 * invArootPi;
270 >      // Half the Smith self piece:
271 >      selfMult1_ = - a2 * invArootPi;
272 >      selfMult2_ = - 2.0 * a2 * a2 * invArootPi / 3.0;
273 >      selfMult4_ = - 4.0 * a2 * a2 * a2 * invArootPi / 5.0;
274      } else {
275        a2 = 0.0;
276        b0c = 1.0 / r;
# Line 271 | Line 279 | namespace OpenMD {
279        b3c = (5.0 * b2c) / r2;
280        b4c = (7.0 * b3c) / r2;
281        b5c = (9.0 * b4c) / r2;
282 <      selfMult_ = b0c;
282 >      selfMult1_ = 0.0;
283 >      selfMult2_ = 0.0;
284 >      selfMult4_ = 0.0;
285      }
286  
287      // higher derivatives of B_0 at the cutoff radius:
# Line 279 | Line 289 | namespace OpenMD {
289      db0c_2 =     -b1c + r2 * b2c;
290      db0c_3 =          3.0*r*b2c  - r2*r*b3c;
291      db0c_4 =          3.0*b2c  - 6.0*r2*b3c     + r2*r2*b4c;
292 <    db0c_5 =                    -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c;
283 <    
292 >    db0c_5 =                    -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c;  
293  
294 +    if (summationMethod_ != esm_EWALD_FULL) {
295 +      selfMult1_ -= b0c;
296 +      selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) /  3.0;
297 +      selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0;
298 +    }
299 +
300      // working variables for the splines:
301      RealType ri, ri2;
302      RealType b0, b1, b2, b3, b4, b5;
# Line 317 | Line 332 | namespace OpenMD {
332      vector<RealType> v31v, v32v;
333      vector<RealType> v41v, v42v, v43v;
334  
320 –    /*
321 –    vector<RealType> dv01v;
322 –    vector<RealType> dv11v;
323 –    vector<RealType> dv21v, dv22v;
324 –    vector<RealType> dv31v, dv32v;
325 –    vector<RealType> dv41v, dv42v, dv43v;
326 –    */
327 –
335      for (int i = 1; i < np_ + 1; i++) {
336        r = RealType(i) * dx;
337        rv.push_back(r);
# Line 488 | Line 495 | namespace OpenMD {
495  
496        case esm_SWITCHING_FUNCTION:
497        case esm_HARD:
498 +      case esm_EWALD_FULL:
499  
500          v01 = f;
501          v11 = g;
# Line 546 | Line 554 | namespace OpenMD {
554  
555          break;
556                  
549 –      case esm_EWALD_FULL:
557        case esm_EWALD_PME:
558        case esm_EWALD_SPME:
559        default :
# Line 575 | Line 582 | namespace OpenMD {
582        v41v.push_back(v41);
583        v42v.push_back(v42);
584        v43v.push_back(v43);
578 –      /*
579 –      dv01v.push_back(dv01);
580 –      dv11v.push_back(dv11);
581 –      dv21v.push_back(dv21);
582 –      dv22v.push_back(dv22);
583 –      dv31v.push_back(dv31);
584 –      dv32v.push_back(dv32);      
585 –      dv41v.push_back(dv41);
586 –      dv42v.push_back(dv42);
587 –      dv43v.push_back(dv43);
588 –      */
585      }
586  
587      // construct the spline structures and fill them with the values we've
# Line 610 | Line 606 | namespace OpenMD {
606      v43s = new CubicSpline();
607      v43s->addPoints(rv, v43v);
608  
613 –    /*
614 –    dv01s = new CubicSpline();
615 –    dv01s->addPoints(rv, dv01v);
616 –    dv11s = new CubicSpline();
617 –    dv11s->addPoints(rv, dv11v);
618 –    dv21s = new CubicSpline();
619 –    dv21s->addPoints(rv, dv21v);
620 –    dv22s = new CubicSpline();
621 –    dv22s->addPoints(rv, dv22v);
622 –    dv31s = new CubicSpline();
623 –    dv31s->addPoints(rv, dv31v);
624 –    dv32s = new CubicSpline();
625 –    dv32s->addPoints(rv, dv32v);
626 –    dv41s = new CubicSpline();
627 –    dv41s->addPoints(rv, dv41v);
628 –    dv42s = new CubicSpline();
629 –    dv42s->addPoints(rv, dv42v);
630 –    dv43s = new CubicSpline();
631 –    dv43s->addPoints(rv, dv43v);
632 –    */
633 –
609      haveElectroSplines_ = true;
610  
611      initialized_ = true;
# Line 704 | Line 679 | namespace OpenMD {
679        FQtids[atid] = fqtid;
680        Jij[fqtid].resize(nFlucq_);
681  
682 <      // Now, iterate over all known fluctuating and add to the coulomb integral map:
682 >      // Now, iterate over all known fluctuating and add to the
683 >      // coulomb integral map:
684        
685        std::set<int>::iterator it;
686        for( it = FQtypes.begin(); it != FQtypes.end(); ++it) {    
# Line 790 | Line 766 | namespace OpenMD {
766      // Excluded potential that is still computed for fluctuating charges
767      excluded_Pot= 0.0;
768  
793 –
769      // some variables we'll need independent of electrostatic type:
770  
771      ri = 1.0 /  *(idat.rij);
# Line 911 | Line 886 | namespace OpenMD {
886          Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or
887                          + rdQbr * rhat * (dv22 - 2.0*v22or));
888      }
889 <    
889 >        
890 >
891      if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
892        J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
893      }    
894 <    
894 >
895      if (a_is_Charge) {    
896        
897        if (b_is_Charge) {
898          pref =  pre11_ * *(idat.electroMult);      
899          U  += C_a * C_b * pref * v01;
900          F  += C_a * C_b * pref * dv01 * rhat;
901 <        
901 >
902          // If this is an excluded pair, there are still indirect
903          // interactions via the reaction field we must worry about:
904  
# Line 931 | Line 907 | namespace OpenMD {
907            indirect_Pot += rfContrib;
908            indirect_F   += rfContrib * 2.0 * ri * rhat;
909          }
910 <        
910 >
911          // Fluctuating charge forces are handled via Coulomb integrals
912          // for excluded pairs (i.e. those connected via bonds) and
913          // with the standard charge-charge interaction otherwise.
914  
915 <        if (idat.excluded) {          
915 >        if (idat.excluded) {
916            if (a_is_Fluctuating || b_is_Fluctuating) {
917              coulInt = J->getValueAt( *(idat.rij) );
918 <            if (a_is_Fluctuating)  dUdCa += coulInt * C_b;
919 <            if (b_is_Fluctuating)  dUdCb += coulInt * C_a;
920 <            excluded_Pot += C_a * C_b * coulInt;
945 <          }          
918 >            if (a_is_Fluctuating) dUdCa += C_b * coulInt;
919 >            if (b_is_Fluctuating) dUdCb += C_a * coulInt;          
920 >          }
921          } else {
922            if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
923 <          if (a_is_Fluctuating) dUdCb += C_a * pref * v01;
924 <        }
923 >          if (b_is_Fluctuating) dUdCb += C_a * pref * v01;
924 >        }              
925        }
926  
927        if (b_is_Dipole) {
# Line 1104 | Line 1079 | namespace OpenMD {
1079  
1080          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1081  
1107 –        //  cerr << " tsum = " << Ta + Tb - cross(  *(idat.d) , F ) << "\n";
1082        }
1083      }
1084  
# Line 1156 | Line 1130 | namespace OpenMD {
1130      // logicals
1131      bool i_is_Charge = data.is_Charge;
1132      bool i_is_Dipole = data.is_Dipole;
1133 +    bool i_is_Quadrupole = data.is_Quadrupole;
1134      bool i_is_Fluctuating = data.is_Fluctuating;
1135      RealType C_a = data.fixedCharge;  
1136 <    RealType self, preVal, DadDa;
1137 <    
1136 >    RealType self(0.0), preVal, DdD, trQ, trQQ;
1137 >
1138 >    if (i_is_Dipole) {
1139 >      DdD = data.dipole.lengthSquare();
1140 >    }
1141 >        
1142      if (i_is_Fluctuating) {
1143        C_a += *(sdat.flucQ);
1144 <      // dVdFQ is really a force, so this is negative the derivative
1166 <      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1144 >      *(sdat.flucQfrc) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1145        (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1146          (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1147      }
# Line 1180 | Line 1158 | namespace OpenMD {
1158        }
1159  
1160        if (i_is_Dipole) {
1161 <        DadDa = data.dipole.lengthSquare();
1184 <        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DadDa;
1161 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DdD;
1162        }
1163        
1164        break;
1165        
1166      case esm_SHIFTED_FORCE:
1167      case esm_SHIFTED_POTENTIAL:
1168 <      if (i_is_Charge) {
1169 <        self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
1170 <        (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1168 >    case esm_TAYLOR_SHIFTED:
1169 >    case esm_EWALD_FULL:
1170 >      if (i_is_Charge)
1171 >        self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge));      
1172 >      if (i_is_Dipole)
1173 >        self += selfMult2_ * pre22_ * DdD;      
1174 >      if (i_is_Quadrupole) {
1175 >        trQ = data.quadrupole.trace();
1176 >        trQQ = (data.quadrupole * data.quadrupole).trace();
1177 >        self += selfMult4_ * pre44_ * (2.0*trQQ + trQ*trQ);
1178 >        if (i_is_Charge)
1179 >          self -= selfMult2_ * pre14_ * 2.0 * C_a * trQ;
1180        }
1181 +      (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;      
1182        break;
1183      default:
1184        break;
# Line 1204 | Line 1191 | namespace OpenMD {
1191      // 12 angstroms seems to be a reasonably good guess for most
1192      // cases.
1193      return 12.0;
1194 +  }
1195 +
1196 +
1197 +  void Electrostatic::ReciprocalSpaceSum(RealType& pot) {
1198 +    
1199 +    RealType kPot = 0.0;
1200 +    RealType kVir = 0.0;
1201 +    
1202 +    const RealType mPoleConverter = 0.20819434; // converts from the
1203 +                                                // internal units of
1204 +                                                // Debye (for dipoles)
1205 +                                                // or Debye-angstroms
1206 +                                                // (for quadrupoles) to
1207 +                                                // electron angstroms or
1208 +                                                // electron-angstroms^2
1209 +    
1210 +    const RealType eConverter = 332.0637778; // convert the
1211 +                                             // Charge-Charge
1212 +                                             // electrostatic
1213 +                                             // interactions into kcal /
1214 +                                             // mol assuming distances
1215 +                                             // are measured in
1216 +                                             // angstroms.
1217 +
1218 +    Mat3x3d hmat = info_->getSnapshotManager()->getCurrentSnapshot()->getHmat();
1219 +    Vector3d box = hmat.diagonals();
1220 +    RealType boxMax = box.max();
1221 +    
1222 +    //int kMax = int(2.0 * M_PI / (pow(dampingAlpha_,2)*cutoffRadius_ * boxMax) );
1223 +    int kMax = 7;
1224 +    int kSqMax = kMax*kMax + 2;
1225 +    
1226 +    int kLimit = kMax+1;
1227 +    int kLim2 = 2*kMax+1;
1228 +    int kSqLim = kSqMax;
1229 +    
1230 +    vector<RealType> AK(kSqLim+1, 0.0);
1231 +    RealType xcl = 2.0 * M_PI / box.x();
1232 +    RealType ycl = 2.0 * M_PI / box.y();
1233 +    RealType zcl = 2.0 * M_PI / box.z();
1234 +    RealType rcl = 2.0 * M_PI / boxMax;
1235 +    RealType rvol = 2.0 * M_PI /(box.x() * box.y() * box.z());
1236 +    
1237 +    if(dampingAlpha_ < 1.0e-12) return;
1238 +    
1239 +    RealType ralph = -0.25/pow(dampingAlpha_,2);
1240 +    
1241 +    // Calculate and store exponential factors  
1242 +    
1243 +    vector<vector<RealType> > elc;
1244 +    vector<vector<RealType> > emc;
1245 +    vector<vector<RealType> > enc;
1246 +    vector<vector<RealType> > els;
1247 +    vector<vector<RealType> > ems;
1248 +    vector<vector<RealType> > ens;
1249 +
1250 +    
1251 +    int nMax = info_->getNAtoms();
1252 +    
1253 +    elc.resize(kLimit+1);
1254 +    emc.resize(kLimit+1);
1255 +    enc.resize(kLimit+1);
1256 +    els.resize(kLimit+1);
1257 +    ems.resize(kLimit+1);
1258 +    ens.resize(kLimit+1);
1259 +
1260 +    for (int j = 0; j < kLimit+1; j++) {
1261 +      elc[j].resize(nMax);
1262 +      emc[j].resize(nMax);
1263 +      enc[j].resize(nMax);
1264 +      els[j].resize(nMax);
1265 +      ems[j].resize(nMax);
1266 +      ens[j].resize(nMax);
1267 +    }
1268 +    
1269 +    Vector3d t( 2.0 * M_PI );
1270 +    t.Vdiv(t, box);
1271 +
1272 +    
1273 +    SimInfo::MoleculeIterator mi;
1274 +    Molecule::AtomIterator ai;
1275 +    int i;
1276 +    Vector3d r;
1277 +    Vector3d tt;
1278 +    
1279 +    for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1280 +         mol = info_->nextMolecule(mi)) {
1281 +      for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1282 +          atom = mol->nextAtom(ai)) {  
1283 +        
1284 +        i = atom->getLocalIndex();
1285 +        r = atom->getPos();
1286 +        info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r);
1287 +        
1288 +        tt.Vmul(t, r);
1289 +
1290 +        elc[1][i] = 1.0;
1291 +        emc[1][i] = 1.0;
1292 +        enc[1][i] = 1.0;
1293 +        els[1][i] = 0.0;
1294 +        ems[1][i] = 0.0;
1295 +        ens[1][i] = 0.0;
1296 +
1297 +        elc[2][i] = cos(tt.x());
1298 +        emc[2][i] = cos(tt.y());
1299 +        enc[2][i] = cos(tt.z());
1300 +        els[2][i] = sin(tt.x());
1301 +        ems[2][i] = sin(tt.y());
1302 +        ens[2][i] = sin(tt.z());
1303 +        
1304 +        for(int l = 3; l <= kLimit; l++) {
1305 +          elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i];
1306 +          emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i];
1307 +          enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i];
1308 +          els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i];
1309 +          ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i];
1310 +          ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i];
1311 +        }
1312 +      }
1313 +    }
1314 +    
1315 +    // Calculate and store AK coefficients:
1316 +    
1317 +    RealType eksq = 1.0;
1318 +    RealType expf = 0.0;
1319 +    if (ralph < 0.0) expf = exp(ralph*rcl*rcl);
1320 +    for (i = 1; i <= kSqLim; i++) {
1321 +      RealType rksq = float(i)*rcl*rcl;
1322 +      eksq = expf*eksq;
1323 +      AK[i] = eConverter * eksq/rksq;
1324 +    }
1325 +    
1326 +    /*
1327 +     * Loop over all k vectors k = 2 pi (ll/Lx, mm/Ly, nn/Lz)
1328 +     * the values of ll, mm and nn are selected so that the symmetry of
1329 +     * reciprocal lattice is taken into account i.e. the following
1330 +     * rules apply.
1331 +     *
1332 +     * ll ranges over the values 0 to kMax only.
1333 +     *
1334 +     * mm ranges over 0 to kMax when ll=0 and over
1335 +     *            -kMax to kMax otherwise.
1336 +     * nn ranges over 1 to kMax when ll=mm=0 and over
1337 +     *            -kMax to kMax otherwise.
1338 +     *
1339 +     * Hence the result of the summation must be doubled at the end.    
1340 +     */
1341 +    
1342 +    std::vector<RealType> clm(nMax, 0.0);
1343 +    std::vector<RealType> slm(nMax, 0.0);
1344 +    std::vector<RealType> ckr(nMax, 0.0);
1345 +    std::vector<RealType> skr(nMax, 0.0);
1346 +    std::vector<RealType> ckc(nMax, 0.0);
1347 +    std::vector<RealType> cks(nMax, 0.0);
1348 +    std::vector<RealType> dkc(nMax, 0.0);
1349 +    std::vector<RealType> dks(nMax, 0.0);
1350 +    std::vector<RealType> qkc(nMax, 0.0);
1351 +    std::vector<RealType> qks(nMax, 0.0);
1352 +    std::vector<Vector3d> dxk(nMax, V3Zero);
1353 +    std::vector<Vector3d> qxk(nMax, V3Zero);
1354 +    RealType rl, rm, rn;
1355 +    Vector3d kVec;
1356 +    Vector3d Qk;
1357 +    Mat3x3d k2;
1358 +    RealType ckcs, ckss, dkcs, dkss, qkcs, qkss;
1359 +    int atid;
1360 +    ElectrostaticAtomData data;
1361 +    RealType C, dk, qk;
1362 +    Vector3d D;
1363 +    Mat3x3d  Q;
1364 +
1365 +    int mMin = kLimit;
1366 +    int nMin = kLimit + 1;
1367 +    for (int l = 1; l <= kLimit; l++) {
1368 +      int ll = l - 1;
1369 +      rl = xcl * float(ll);
1370 +      for (int mmm = mMin; mmm <= kLim2; mmm++) {
1371 +        int mm = mmm - kLimit;
1372 +        int m = abs(mm) + 1;
1373 +        rm = ycl * float(mm);
1374 +        // Set temporary products of exponential terms
1375 +        for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1376 +             mol = info_->nextMolecule(mi)) {
1377 +          for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1378 +              atom = mol->nextAtom(ai)) {
1379 +            
1380 +            i = atom->getLocalIndex();
1381 +            if(mm < 0) {
1382 +              clm[i]=elc[l][i]*emc[m][i]+els[l][i]*ems[m][i];
1383 +              slm[i]=els[l][i]*emc[m][i]-ems[m][i]*elc[l][i];
1384 +            } else {
1385 +              clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i];
1386 +              slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i];
1387 +            }
1388 +          }
1389 +        }
1390 +        for (int nnn = nMin; nnn <= kLim2; nnn++) {
1391 +          int nn = nnn - kLimit;          
1392 +          int n = abs(nn) + 1;
1393 +          rn = zcl * float(nn);
1394 +          // Test on magnitude of k vector:
1395 +          int kk=ll*ll + mm*mm + nn*nn;
1396 +          if(kk <= kSqLim) {
1397 +            kVec = Vector3d(rl, rm, rn);
1398 +            k2 = outProduct(kVec, kVec);
1399 +            // Calculate exp(ikr) terms
1400 +            for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1401 +                 mol = info_->nextMolecule(mi)) {
1402 +              for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1403 +                  atom = mol->nextAtom(ai)) {
1404 +                i = atom->getLocalIndex();
1405 +                
1406 +                if (nn < 0) {
1407 +                  ckr[i]=clm[i]*enc[n][i]+slm[i]*ens[n][i];
1408 +                  skr[i]=slm[i]*enc[n][i]-clm[i]*ens[n][i];
1409 +
1410 +                } else {
1411 +                  ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i];
1412 +                  skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i];
1413 +                }
1414 +              }
1415 +            }
1416 +            
1417 +            // Calculate scalar and vector products for each site:
1418 +            
1419 +            for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1420 +                 mol = info_->nextMolecule(mi)) {
1421 +              for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1422 +                  atom = mol->nextAtom(ai)) {
1423 +                i = atom->getLocalIndex();
1424 +                int atid = atom->getAtomType()->getIdent();
1425 +                data = ElectrostaticMap[Etids[atid]];
1426 +                              
1427 +                if (data.is_Charge) {
1428 +                  C = data.fixedCharge;
1429 +                  if (atom->isFluctuatingCharge()) C += atom->getFlucQPos();
1430 +                  ckc[i] = C * ckr[i];
1431 +                  cks[i] = C * skr[i];
1432 +                }
1433 +                
1434 +                if (data.is_Dipole) {
1435 +                  D = atom->getDipole() * mPoleConverter;
1436 +                  dk = dot(D, kVec);
1437 +                  dxk[i] = cross(D, kVec);
1438 +                  dkc[i] = dk * ckr[i];
1439 +                  dks[i] = dk * skr[i];
1440 +                }
1441 +                if (data.is_Quadrupole) {
1442 +                  Q = atom->getQuadrupole() * mPoleConverter;
1443 +                  Qk = Q * kVec;                  
1444 +                  qk = dot(Qk, kVec);
1445 +                  qxk[i] = cross(Qk, kVec);
1446 +                  qkc[i] = qk * ckr[i];
1447 +                  qks[i] = qk * skr[i];
1448 +                }              
1449 +              }
1450 +            }
1451 +
1452 +            // calculate vector sums
1453 +            
1454 +            ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1455 +            ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1456 +            dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1457 +            dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1458 +            qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1459 +            qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1460 +            
1461 + #ifdef IS_MPI
1462 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE,
1463 +                                      MPI::SUM);
1464 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckss, 1, MPI::REALTYPE,
1465 +                                      MPI::SUM);
1466 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkcs, 1, MPI::REALTYPE,
1467 +                                      MPI::SUM);
1468 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkss, 1, MPI::REALTYPE,
1469 +                                      MPI::SUM);
1470 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkcs, 1, MPI::REALTYPE,
1471 +                                      MPI::SUM);
1472 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkss, 1, MPI::REALTYPE,
1473 +                                      MPI::SUM);
1474 + #endif        
1475 +            
1476 +            // Accumulate potential energy and virial contribution:
1477 +
1478 +            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1479 +                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs));
1480 +
1481 +            kVir += 2.0 * rvol  * AK[kk]*(ckcs*ckcs+ckss*ckss
1482 +                                          +4.0*(ckss*dkcs-ckcs*dkss)
1483 +                                          +3.0*(dkcs*dkcs+dkss*dkss)
1484 +                                          -6.0*(ckss*qkss+ckcs*qkcs)
1485 +                                          +8.0*(dkss*qkcs-dkcs*qkss)
1486 +                                          +5.0*(qkss*qkss+qkcs*qkcs));
1487 +            
1488 +            // Calculate force and torque for each site:
1489 +            
1490 +            for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1491 +                 mol = info_->nextMolecule(mi)) {
1492 +              for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1493 +                  atom = mol->nextAtom(ai)) {
1494 +                
1495 +                i = atom->getLocalIndex();
1496 +                atid = atom->getAtomType()->getIdent();
1497 +                data = ElectrostaticMap[Etids[atid]];
1498 +
1499 +                RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs)
1500 +                                     - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss));
1501 +                RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
1502 +                                         -ckr[i]*(ckss+dkcs-qkss));
1503 +                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
1504 +                                             +skr[i]*(ckss+dkcs-qkss));
1505 +              
1506 +                atom->addFrc( 4.0 * rvol * qfrc * kVec );
1507 +                
1508 +                if (data.is_Dipole) {
1509 +                  atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] );
1510 +                }
1511 +                if (data.is_Quadrupole) {
1512 +                  atom->addTrq( 4.0 * rvol * qtrq2 * qxk[i] );
1513 +                }
1514 +              }
1515 +            }
1516 +          }
1517 +        }
1518 +        nMin = 1;
1519 +      }
1520 +      mMin = 1;
1521 +    }
1522 +    pot += kPot;  
1523    }
1524   }

Diff Legend

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