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 1907 by gezelter, Fri Jul 19 18:18:27 2013 UTC vs.
Revision 1993 by gezelter, Tue Apr 29 17:32:31 2014 UTC

# Line 40 | Line 40
40   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43 + #ifdef IS_MPI
44 + #include <mpi.h>
45 + #endif
46 +
47   #include <stdio.h>
48   #include <string.h>
49  
50   #include <cmath>
51 + #include <numeric>
52   #include "nonbonded/Electrostatic.hpp"
53   #include "utils/simError.h"
54   #include "types/NonBondedInteractionType.hpp"
# Line 55 | Line 60
60   #include "utils/PhysicalConstants.hpp"
61   #include "math/erfc.hpp"
62   #include "math/SquareMatrix.hpp"
63 + #include "primitives/Molecule.hpp"
64 + #include "flucq/FluctuatingChargeForces.hpp"
65  
66   namespace OpenMD {
67    
# Line 64 | Line 71 | namespace OpenMD {
71                                    haveDampingAlpha_(false),
72                                    haveDielectric_(false),
73                                    haveElectroSplines_(false)
74 <  {}
74 >  {
75 >    flucQ_ = new FluctuatingChargeForces(info_);
76 >  }
77    
78 +  void Electrostatic::setForceField(ForceField *ff) {
79 +    forceField_ = ff;
80 +    flucQ_->setForceField(forceField_);
81 +  }
82 +
83 +  void Electrostatic::setSimulatedAtomTypes(set<AtomType*> &simtypes) {
84 +    simTypes_ = simtypes;
85 +    flucQ_->setSimulatedAtomTypes(simTypes_);
86 +  }
87 +
88    void Electrostatic::initialize() {
89      
90      Globals* simParams_ = info_->getSimParams();
# Line 191 | Line 210 | namespace OpenMD {
210        simError();
211      }
212            
213 <    if (screeningMethod_ == DAMPED) {      
213 >    if (screeningMethod_ == DAMPED || summationMethod_ == esm_EWALD_FULL) {
214        if (!simParams_->haveDampingAlpha()) {
215          // first set a cutoff dependent alpha value
216          // we assume alpha depends linearly with rcut from 0 to 20.5 ang
217          dampingAlpha_ = 0.425 - cutoffRadius_* 0.02;
218 <        if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;
200 <        
218 >        if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;        
219          // throw warning
220          sprintf( painCave.errMsg,
221                   "Electrostatic::initialize: dampingAlpha was not specified in the\n"
# Line 213 | Line 231 | namespace OpenMD {
231        haveDampingAlpha_ = true;
232      }
233  
234 +
235      Etypes.clear();
236      Etids.clear();
237      FQtypes.clear();
# Line 262 | Line 281 | namespace OpenMD {
281        b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2;
282        b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2;
283        b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2;
265 –      //selfMult1_ = - 2.0 * a2 * invArootPi;
266 –      //selfMult2_ = - 4.0 * a2 * a2 * invArootPi / 3.0;
267 –      //selfMult4_ = - 8.0 * a2 * a2 * a2 * invArootPi / 5.0;
284        // Half the Smith self piece:
285        selfMult1_ = - a2 * invArootPi;
286        selfMult2_ = - 2.0 * a2 * a2 * invArootPi / 3.0;
# Line 288 | Line 304 | namespace OpenMD {
304      db0c_3 =          3.0*r*b2c  - r2*r*b3c;
305      db0c_4 =          3.0*b2c  - 6.0*r2*b3c     + r2*r2*b4c;
306      db0c_5 =                    -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c;  
291 –    
292 –    selfMult1_ -= b0c;
293 –    selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) /  3.0;
294 –    selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0;
307  
308 +    if (summationMethod_ != esm_EWALD_FULL) {
309 +      selfMult1_ -= b0c;
310 +      selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) /  3.0;
311 +      selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0;
312 +    }
313 +
314      // working variables for the splines:
315      RealType ri, ri2;
316      RealType b0, b1, b2, b3, b4, b5;
# Line 328 | Line 346 | namespace OpenMD {
346      vector<RealType> v31v, v32v;
347      vector<RealType> v41v, v42v, v43v;
348  
331 –    /*
332 –    vector<RealType> dv01v;
333 –    vector<RealType> dv11v;
334 –    vector<RealType> dv21v, dv22v;
335 –    vector<RealType> dv31v, dv32v;
336 –    vector<RealType> dv41v, dv42v, dv43v;
337 –    */
338 –
349      for (int i = 1; i < np_ + 1; i++) {
350        r = RealType(i) * dx;
351        rv.push_back(r);
# Line 499 | Line 509 | namespace OpenMD {
509  
510        case esm_SWITCHING_FUNCTION:
511        case esm_HARD:
512 +      case esm_EWALD_FULL:
513  
514          v01 = f;
515          v11 = g;
# Line 557 | Line 568 | namespace OpenMD {
568  
569          break;
570                  
560 –      case esm_EWALD_FULL:
571        case esm_EWALD_PME:
572        case esm_EWALD_SPME:
573        default :
# Line 586 | Line 596 | namespace OpenMD {
596        v41v.push_back(v41);
597        v42v.push_back(v42);
598        v43v.push_back(v43);
589 –      /*
590 –      dv01v.push_back(dv01);
591 –      dv11v.push_back(dv11);
592 –      dv21v.push_back(dv21);
593 –      dv22v.push_back(dv22);
594 –      dv31v.push_back(dv31);
595 –      dv32v.push_back(dv32);      
596 –      dv41v.push_back(dv41);
597 –      dv42v.push_back(dv42);
598 –      dv43v.push_back(dv43);
599 –      */
599      }
600  
601      // construct the spline structures and fill them with the values we've
# Line 621 | Line 620 | namespace OpenMD {
620      v43s = new CubicSpline();
621      v43s->addPoints(rv, v43v);
622  
624 –    /*
625 –    dv01s = new CubicSpline();
626 –    dv01s->addPoints(rv, dv01v);
627 –    dv11s = new CubicSpline();
628 –    dv11s->addPoints(rv, dv11v);
629 –    dv21s = new CubicSpline();
630 –    dv21s->addPoints(rv, dv21v);
631 –    dv22s = new CubicSpline();
632 –    dv22s->addPoints(rv, dv22v);
633 –    dv31s = new CubicSpline();
634 –    dv31s->addPoints(rv, dv31v);
635 –    dv32s = new CubicSpline();
636 –    dv32s->addPoints(rv, dv32v);
637 –    dv41s = new CubicSpline();
638 –    dv41s->addPoints(rv, dv41v);
639 –    dv42s = new CubicSpline();
640 –    dv42s->addPoints(rv, dv42v);
641 –    dv43s = new CubicSpline();
642 –    dv43s->addPoints(rv, dv43v);
643 –    */
644 –
623      haveElectroSplines_ = true;
624  
625      initialized_ = true;
# Line 715 | Line 693 | namespace OpenMD {
693        FQtids[atid] = fqtid;
694        Jij[fqtid].resize(nFlucq_);
695  
696 <      // Now, iterate over all known fluctuating and add to the coulomb integral map:
696 >      // Now, iterate over all known fluctuating and add to the
697 >      // coulomb integral map:
698        
699        std::set<int>::iterator it;
700        for( it = FQtypes.begin(); it != FQtypes.end(); ++it) {    
# Line 789 | Line 768 | namespace OpenMD {
768      Tb.zero(); // Torque on site b
769      Ea.zero(); // Electric field at site a
770      Eb.zero(); // Electric field at site b
771 +    Pa = 0.0;  // Site potential at site a
772 +    Pb = 0.0;  // Site potential at site b
773      dUdCa = 0.0; // fluctuating charge force at site a
774      dUdCb = 0.0; // fluctuating charge force at site a
775      
# Line 801 | Line 782 | namespace OpenMD {
782      // Excluded potential that is still computed for fluctuating charges
783      excluded_Pot= 0.0;
784  
804 –
785      // some variables we'll need independent of electrostatic type:
786  
787      ri = 1.0 /  *(idat.rij);
# Line 864 | Line 844 | namespace OpenMD {
844        if (idat.excluded) {
845          *(idat.skippedCharge2) += C_a;
846        } else {
847 <        // only do the field if we're not excluded:
847 >        // only do the field and site potentials if we're not excluded:
848          Eb -= C_a *  pre11_ * dv01 * rhat;
849 +        Pb += C_a *  pre11_ * v01;
850        }
851      }
852      
# Line 873 | Line 854 | namespace OpenMD {
854        D_a = *(idat.dipole1);
855        rdDa = dot(rhat, D_a);
856        rxDa = cross(rhat, D_a);
857 <      if (!idat.excluded)
857 >      if (!idat.excluded) {
858          Eb -=  pre12_ * ((dv11-v11or) * rdDa * rhat + v11or * D_a);
859 +        Pb +=  pre12_ * v11 * rdDa;
860 +      }
861 +
862      }
863      
864      if (a_is_Quadrupole) {
# Line 884 | Line 868 | namespace OpenMD {
868        rQa = rhat * Q_a;
869        rdQar = dot(rhat, Qar);
870        rxQar = cross(rhat, Qar);
871 <      if (!idat.excluded)
871 >      if (!idat.excluded) {
872          Eb -= pre14_ * (trQa * rhat * dv21 + 2.0 * Qar * v22or
873                          + rdQar * rhat * (dv22 - 2.0*v22or));
874 +        Pb += pre14_ * (v21 * trQa + v22 * rdQar);
875 +      }
876      }
877      
878      if (b_is_Charge) {
# Line 900 | Line 886 | namespace OpenMD {
886        } else {
887          // only do the field if we're not excluded:
888          Ea += C_b *  pre11_ * dv01 * rhat;
889 +        Pa += C_b *  pre11_ * v01;
890 +
891        }
892      }
893      
# Line 907 | Line 895 | namespace OpenMD {
895        D_b = *(idat.dipole2);
896        rdDb = dot(rhat, D_b);
897        rxDb = cross(rhat, D_b);
898 <      if (!idat.excluded)
898 >      if (!idat.excluded) {
899          Ea += pre12_ * ((dv11-v11or) * rdDb * rhat + v11or * D_b);
900 +        Pa += pre12_ * v11 * rdDb;
901 +      }
902      }
903      
904      if (b_is_Quadrupole) {
# Line 918 | Line 908 | namespace OpenMD {
908        rQb = rhat * Q_b;
909        rdQbr = dot(rhat, Qbr);
910        rxQbr = cross(rhat, Qbr);
911 <      if (!idat.excluded)
911 >      if (!idat.excluded) {
912          Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or
913                          + rdQbr * rhat * (dv22 - 2.0*v22or));
914 +        Pa += pre14_ * (v21 * trQb + v22 * rdQbr);
915 +      }
916      }
917 <    
917 >        
918 >
919      if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
920        J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
921      }    
922 <    
922 >
923      if (a_is_Charge) {    
924        
925        if (b_is_Charge) {
926          pref =  pre11_ * *(idat.electroMult);      
927          U  += C_a * C_b * pref * v01;
928          F  += C_a * C_b * pref * dv01 * rhat;
929 <        
929 >
930          // If this is an excluded pair, there are still indirect
931          // interactions via the reaction field we must worry about:
932  
# Line 942 | Line 935 | namespace OpenMD {
935            indirect_Pot += rfContrib;
936            indirect_F   += rfContrib * 2.0 * ri * rhat;
937          }
938 <        
938 >
939          // Fluctuating charge forces are handled via Coulomb integrals
940          // for excluded pairs (i.e. those connected via bonds) and
941          // with the standard charge-charge interaction otherwise.
942  
943 <        if (idat.excluded) {          
943 >        if (idat.excluded) {
944            if (a_is_Fluctuating || b_is_Fluctuating) {
945              coulInt = J->getValueAt( *(idat.rij) );
946 <            if (a_is_Fluctuating)  dUdCa += coulInt * C_b;
947 <            if (b_is_Fluctuating)  dUdCb += coulInt * C_a;
948 <            excluded_Pot += C_a * C_b * coulInt;
956 <          }          
946 >            if (a_is_Fluctuating) dUdCa += C_b * coulInt;
947 >            if (b_is_Fluctuating) dUdCb += C_a * coulInt;          
948 >          }
949          } else {
950            if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
951 <          if (a_is_Fluctuating) dUdCb += C_a * pref * v01;
952 <        }
951 >          if (b_is_Fluctuating) dUdCb += C_a * pref * v01;
952 >        }              
953        }
954  
955        if (b_is_Dipole) {
# Line 1023 | Line 1015 | namespace OpenMD {
1015          F  -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat;
1016          Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
1017          Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
1026 –
1018          // Even if we excluded this pair from direct interactions, we
1019          // still have the reaction-field-mediated dipole-dipole
1020          // interaction:
# Line 1083 | Line 1074 | namespace OpenMD {
1074          trQaQb = QaQb.trace();
1075          rQaQb = rhat * QaQb;
1076          QaQbr = QaQb * rhat;
1077 <        QaxQb = cross(Q_a, Q_b);
1077 >        QaxQb = mCross(Q_a, Q_b);
1078          rQaQbr = dot(rQa, Qbr);
1079          rQaxQbr = cross(rQa, Qbr);
1080          
# Line 1114 | Line 1105 | namespace OpenMD {
1105          //             + 4.0 * cross(rhat, QbQar)
1106  
1107          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1117 –
1108        }
1109      }
1110  
# Line 1123 | Line 1113 | namespace OpenMD {
1113        *(idat.eField2) += Eb * *(idat.electroMult);
1114      }
1115  
1116 +    if (idat.doSitePotential) {
1117 +      *(idat.sPot1) += Pa * *(idat.electroMult);
1118 +      *(idat.sPot2) += Pb * *(idat.electroMult);
1119 +    }
1120 +
1121      if (a_is_Fluctuating) *(idat.dVdFQ1) += dUdCa * *(idat.sw);
1122      if (b_is_Fluctuating) *(idat.dVdFQ2) += dUdCb * *(idat.sw);
1123  
# Line 1177 | Line 1172 | namespace OpenMD {
1172          
1173      if (i_is_Fluctuating) {
1174        C_a += *(sdat.flucQ);
1175 <      // dVdFQ is really a force, so this is negative the derivative
1176 <      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1177 <      (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1178 <        (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1175 >
1176 >      flucQ_->getSelfInteraction(sdat.atid, *(sdat.flucQ),  
1177 >                                 (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY],
1178 >                                 *(sdat.flucQfrc) );
1179 >
1180      }
1181  
1182      switch (summationMethod_) {
# Line 1203 | Line 1199 | namespace OpenMD {
1199      case esm_SHIFTED_FORCE:
1200      case esm_SHIFTED_POTENTIAL:
1201      case esm_TAYLOR_SHIFTED:
1202 +    case esm_EWALD_FULL:
1203        if (i_is_Charge)
1204          self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge));      
1205        if (i_is_Dipole)
# Line 1227 | Line 1224 | namespace OpenMD {
1224      // 12 angstroms seems to be a reasonably good guess for most
1225      // cases.
1226      return 12.0;
1227 +  }
1228 +
1229 +
1230 +  void Electrostatic::ReciprocalSpaceSum(RealType& pot) {
1231 +    
1232 +    RealType kPot = 0.0;
1233 +    RealType kVir = 0.0;
1234 +    
1235 +    const RealType mPoleConverter = 0.20819434; // converts from the
1236 +                                                // internal units of
1237 +                                                // Debye (for dipoles)
1238 +                                                // or Debye-angstroms
1239 +                                                // (for quadrupoles) to
1240 +                                                // electron angstroms or
1241 +                                                // electron-angstroms^2
1242 +    
1243 +    const RealType eConverter = 332.0637778; // convert the
1244 +                                             // Charge-Charge
1245 +                                             // electrostatic
1246 +                                             // interactions into kcal /
1247 +                                             // mol assuming distances
1248 +                                             // are measured in
1249 +                                             // angstroms.
1250 +
1251 +    Mat3x3d hmat = info_->getSnapshotManager()->getCurrentSnapshot()->getHmat();
1252 +    Vector3d box = hmat.diagonals();
1253 +    RealType boxMax = box.max();
1254 +    
1255 +    //int kMax = int(2.0 * M_PI / (pow(dampingAlpha_,2)*cutoffRadius_ * boxMax) );
1256 +    int kMax = 7;
1257 +    int kSqMax = kMax*kMax + 2;
1258 +    
1259 +    int kLimit = kMax+1;
1260 +    int kLim2 = 2*kMax+1;
1261 +    int kSqLim = kSqMax;
1262 +    
1263 +    vector<RealType> AK(kSqLim+1, 0.0);
1264 +    RealType xcl = 2.0 * M_PI / box.x();
1265 +    RealType ycl = 2.0 * M_PI / box.y();
1266 +    RealType zcl = 2.0 * M_PI / box.z();
1267 +    RealType rcl = 2.0 * M_PI / boxMax;
1268 +    RealType rvol = 2.0 * M_PI /(box.x() * box.y() * box.z());
1269 +    
1270 +    if(dampingAlpha_ < 1.0e-12) return;
1271 +    
1272 +    RealType ralph = -0.25/pow(dampingAlpha_,2);
1273 +    
1274 +    // Calculate and store exponential factors  
1275 +    
1276 +    vector<vector<RealType> > elc;
1277 +    vector<vector<RealType> > emc;
1278 +    vector<vector<RealType> > enc;
1279 +    vector<vector<RealType> > els;
1280 +    vector<vector<RealType> > ems;
1281 +    vector<vector<RealType> > ens;
1282 +    
1283 +    int nMax = info_->getNAtoms();
1284 +    
1285 +    elc.resize(kLimit+1);
1286 +    emc.resize(kLimit+1);
1287 +    enc.resize(kLimit+1);
1288 +    els.resize(kLimit+1);
1289 +    ems.resize(kLimit+1);
1290 +    ens.resize(kLimit+1);
1291 +
1292 +    for (int j = 0; j < kLimit+1; j++) {
1293 +      elc[j].resize(nMax);
1294 +      emc[j].resize(nMax);
1295 +      enc[j].resize(nMax);
1296 +      els[j].resize(nMax);
1297 +      ems[j].resize(nMax);
1298 +      ens[j].resize(nMax);
1299 +    }
1300 +    
1301 +    Vector3d t( 2.0 * M_PI );
1302 +    t.Vdiv(t, box);
1303 +
1304 +    SimInfo::MoleculeIterator mi;
1305 +    Molecule::AtomIterator ai;
1306 +    int i;
1307 +    Vector3d r;
1308 +    Vector3d tt;
1309 +    
1310 +    for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1311 +         mol = info_->nextMolecule(mi)) {
1312 +      for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1313 +          atom = mol->nextAtom(ai)) {  
1314 +        
1315 +        i = atom->getLocalIndex();
1316 +        r = atom->getPos();
1317 +        info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r);
1318 +        
1319 +        tt.Vmul(t, r);
1320 +
1321 +        elc[1][i] = 1.0;
1322 +        emc[1][i] = 1.0;
1323 +        enc[1][i] = 1.0;
1324 +        els[1][i] = 0.0;
1325 +        ems[1][i] = 0.0;
1326 +        ens[1][i] = 0.0;
1327 +
1328 +        elc[2][i] = cos(tt.x());
1329 +        emc[2][i] = cos(tt.y());
1330 +        enc[2][i] = cos(tt.z());
1331 +        els[2][i] = sin(tt.x());
1332 +        ems[2][i] = sin(tt.y());
1333 +        ens[2][i] = sin(tt.z());
1334 +        
1335 +        for(int l = 3; l <= kLimit; l++) {
1336 +          elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i];
1337 +          emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i];
1338 +          enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i];
1339 +          els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i];
1340 +          ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i];
1341 +          ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i];
1342 +        }
1343 +      }
1344 +    }
1345 +    
1346 +    // Calculate and store AK coefficients:
1347 +    
1348 +    RealType eksq = 1.0;
1349 +    RealType expf = 0.0;
1350 +    if (ralph < 0.0) expf = exp(ralph*rcl*rcl);
1351 +    for (i = 1; i <= kSqLim; i++) {
1352 +      RealType rksq = float(i)*rcl*rcl;
1353 +      eksq = expf*eksq;
1354 +      AK[i] = eConverter * eksq/rksq;
1355 +    }
1356 +    
1357 +    /*
1358 +     * Loop over all k vectors k = 2 pi (ll/Lx, mm/Ly, nn/Lz)
1359 +     * the values of ll, mm and nn are selected so that the symmetry of
1360 +     * reciprocal lattice is taken into account i.e. the following
1361 +     * rules apply.
1362 +     *
1363 +     * ll ranges over the values 0 to kMax only.
1364 +     *
1365 +     * mm ranges over 0 to kMax when ll=0 and over
1366 +     *            -kMax to kMax otherwise.
1367 +     * nn ranges over 1 to kMax when ll=mm=0 and over
1368 +     *            -kMax to kMax otherwise.
1369 +     *
1370 +     * Hence the result of the summation must be doubled at the end.    
1371 +     */
1372 +    
1373 +    std::vector<RealType> clm(nMax, 0.0);
1374 +    std::vector<RealType> slm(nMax, 0.0);
1375 +    std::vector<RealType> ckr(nMax, 0.0);
1376 +    std::vector<RealType> skr(nMax, 0.0);
1377 +    std::vector<RealType> ckc(nMax, 0.0);
1378 +    std::vector<RealType> cks(nMax, 0.0);
1379 +    std::vector<RealType> dkc(nMax, 0.0);
1380 +    std::vector<RealType> dks(nMax, 0.0);
1381 +    std::vector<RealType> qkc(nMax, 0.0);
1382 +    std::vector<RealType> qks(nMax, 0.0);
1383 +    std::vector<Vector3d> dxk(nMax, V3Zero);
1384 +    std::vector<Vector3d> qxk(nMax, V3Zero);
1385 +    RealType rl, rm, rn;
1386 +    Vector3d kVec;
1387 +    Vector3d Qk;
1388 +    Mat3x3d k2;
1389 +    RealType ckcs, ckss, dkcs, dkss, qkcs, qkss;
1390 +    int atid;
1391 +    ElectrostaticAtomData data;
1392 +    RealType C, dk, qk;
1393 +    Vector3d D;
1394 +    Mat3x3d  Q;
1395 +
1396 +    int mMin = kLimit;
1397 +    int nMin = kLimit + 1;
1398 +    for (int l = 1; l <= kLimit; l++) {
1399 +      int ll = l - 1;
1400 +      rl = xcl * float(ll);
1401 +      for (int mmm = mMin; mmm <= kLim2; mmm++) {
1402 +        int mm = mmm - kLimit;
1403 +        int m = abs(mm) + 1;
1404 +        rm = ycl * float(mm);
1405 +        // Set temporary products of exponential terms
1406 +        for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1407 +             mol = info_->nextMolecule(mi)) {
1408 +          for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1409 +              atom = mol->nextAtom(ai)) {
1410 +            
1411 +            i = atom->getLocalIndex();
1412 +            if(mm < 0) {
1413 +              clm[i]=elc[l][i]*emc[m][i]+els[l][i]*ems[m][i];
1414 +              slm[i]=els[l][i]*emc[m][i]-ems[m][i]*elc[l][i];
1415 +            } else {
1416 +              clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i];
1417 +              slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i];
1418 +            }
1419 +          }
1420 +        }
1421 +        for (int nnn = nMin; nnn <= kLim2; nnn++) {
1422 +          int nn = nnn - kLimit;          
1423 +          int n = abs(nn) + 1;
1424 +          rn = zcl * float(nn);
1425 +          // Test on magnitude of k vector:
1426 +          int kk=ll*ll + mm*mm + nn*nn;
1427 +          if(kk <= kSqLim) {
1428 +            kVec = Vector3d(rl, rm, rn);
1429 +            k2 = outProduct(kVec, kVec);
1430 +            // Calculate exp(ikr) terms
1431 +            for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1432 +                 mol = info_->nextMolecule(mi)) {
1433 +              for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1434 +                  atom = mol->nextAtom(ai)) {
1435 +                i = atom->getLocalIndex();
1436 +                
1437 +                if (nn < 0) {
1438 +                  ckr[i]=clm[i]*enc[n][i]+slm[i]*ens[n][i];
1439 +                  skr[i]=slm[i]*enc[n][i]-clm[i]*ens[n][i];
1440 +
1441 +                } else {
1442 +                  ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i];
1443 +                  skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i];
1444 +                }
1445 +              }
1446 +            }
1447 +            
1448 +            // Calculate scalar and vector products for each site:
1449 +            
1450 +            for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1451 +                 mol = info_->nextMolecule(mi)) {
1452 +              for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1453 +                  atom = mol->nextAtom(ai)) {
1454 +                i = atom->getLocalIndex();
1455 +                int atid = atom->getAtomType()->getIdent();
1456 +                data = ElectrostaticMap[Etids[atid]];
1457 +                              
1458 +                if (data.is_Charge) {
1459 +                  C = data.fixedCharge;
1460 +                  if (atom->isFluctuatingCharge()) C += atom->getFlucQPos();
1461 +                  ckc[i] = C * ckr[i];
1462 +                  cks[i] = C * skr[i];
1463 +                }
1464 +                
1465 +                if (data.is_Dipole) {
1466 +                  D = atom->getDipole() * mPoleConverter;
1467 +                  dk = dot(D, kVec);
1468 +                  dxk[i] = cross(D, kVec);
1469 +                  dkc[i] = dk * ckr[i];
1470 +                  dks[i] = dk * skr[i];
1471 +                }
1472 +                if (data.is_Quadrupole) {
1473 +                  Q = atom->getQuadrupole() * mPoleConverter;
1474 +                  Qk = Q * kVec;                  
1475 +                  qk = dot(kVec, Qk);
1476 +                  qxk[i] = -cross(kVec, Qk);
1477 +                  qkc[i] = qk * ckr[i];
1478 +                  qks[i] = qk * skr[i];
1479 +                }              
1480 +              }
1481 +            }
1482 +
1483 +            // calculate vector sums
1484 +            
1485 +            ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1486 +            ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1487 +            dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1488 +            dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1489 +            qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1490 +            qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1491 +            
1492 + #ifdef IS_MPI
1493 +            MPI_Allreduce(MPI_IN_PLACE, &ckcs, 1, MPI_REALTYPE,
1494 +                          MPI_SUM, MPI_COMM_WORLD);
1495 +            MPI_Allreduce(MPI_IN_PLACE, &ckss, 1, MPI_REALTYPE,
1496 +                          MPI_SUM, MPI_COMM_WORLD);
1497 +            MPI_Allreduce(MPI_IN_PLACE, &dkcs, 1, MPI_REALTYPE,
1498 +                          MPI_SUM, MPI_COMM_WORLD);
1499 +            MPI_Allreduce(MPI_IN_PLACE, &dkss, 1, MPI_REALTYPE,
1500 +                          MPI_SUM, MPI_COMM_WORLD);
1501 +            MPI_Allreduce(MPI_IN_PLACE, &qkcs, 1, MPI_REALTYPE,
1502 +                          MPI_SUM, MPI_COMM_WORLD);
1503 +            MPI_Allreduce(MPI_IN_PLACE, &qkss, 1, MPI_REALTYPE,
1504 +                          MPI_SUM, MPI_COMM_WORLD);
1505 + #endif        
1506 +            
1507 +            // Accumulate potential energy and virial contribution:
1508 +
1509 +            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1510 +                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs));
1511 +
1512 +            kVir += 2.0 * rvol  * AK[kk]*(ckcs*ckcs+ckss*ckss
1513 +                                          +4.0*(ckss*dkcs-ckcs*dkss)
1514 +                                          +3.0*(dkcs*dkcs+dkss*dkss)
1515 +                                          -6.0*(ckss*qkss+ckcs*qkcs)
1516 +                                          +8.0*(dkss*qkcs-dkcs*qkss)
1517 +                                          +5.0*(qkss*qkss+qkcs*qkcs));
1518 +            
1519 +            // Calculate force and torque for each site:
1520 +            
1521 +            for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1522 +                 mol = info_->nextMolecule(mi)) {
1523 +              for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1524 +                  atom = mol->nextAtom(ai)) {
1525 +                
1526 +                i = atom->getLocalIndex();
1527 +                atid = atom->getAtomType()->getIdent();
1528 +                data = ElectrostaticMap[Etids[atid]];
1529 +
1530 +                RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs)
1531 +                                     - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss));
1532 +                RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
1533 +                                         -ckr[i]*(ckss+dkcs-qkss));
1534 +                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
1535 +                                            +skr[i]*(ckss+dkcs-qkss));
1536 +              
1537 +                atom->addFrc( 4.0 * rvol * qfrc * kVec );
1538 +
1539 +                if (atom->isFluctuatingCharge()) {
1540 +                  atom->addFlucQFrc( - 2.0 * rvol * qtrq2 );
1541 +                }
1542 +                  
1543 +                if (data.is_Dipole) {
1544 +                  atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] );
1545 +                }
1546 +                if (data.is_Quadrupole) {
1547 +                  atom->addTrq( 4.0 * rvol * qtrq2 * qxk[i] );
1548 +                }
1549 +              }
1550 +            }
1551 +          }
1552 +        }
1553 +        nMin = 1;
1554 +      }
1555 +      mMin = 1;
1556 +    }
1557 +    pot += kPot;  
1558    }
1559   }

Diff Legend

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