ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/devel_omp/src/brains/ForceManager.cpp
Revision: 1608
Committed: Tue Aug 9 01:58:56 2011 UTC (13 years, 8 months ago) by mciznick
File size: 42806 byte(s)
Log Message: 
First OpenMP version.

File Contents

# Content
1 /*
2 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 *
4 * The University of Notre Dame grants you ("Licensee") a
5 * non-exclusive, royalty free, license to use, modify and
6 * redistribute this software in source and binary code form, provided
7 * that the following conditions are met:
8 *
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 *
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the
15 * distribution.
16 *
17 * This software is provided "AS IS," without a warranty of any
18 * kind. All express or implied conditions, representations and
19 * warranties, including any implied warranty of merchantability,
20 * fitness for a particular purpose or non-infringement, are hereby
21 * excluded. The University of Notre Dame and its licensors shall not
22 * be liable for any damages suffered by licensee as a result of
23 * using, modifying or distributing the software or its
24 * derivatives. In no event will the University of Notre Dame or its
25 * licensors be liable for any lost revenue, profit or data, or for
26 * direct, indirect, special, consequential, incidental or punitive
27 * damages, however caused and regardless of the theory of liability,
28 * arising out of the use of or inability to use software, even if the
29 * University of Notre Dame has been advised of the possibility of
30 * such damages.
31 *
32 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33 * research, please cite the appropriate papers when you publish your
34 * work. Good starting points are:
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).
39 * [4] Vardeman & Gezelter, in progress (2009).
40 */
41
42 /**
43 * @file ForceManager.cpp
44 * @author tlin
45 * @date 11/09/2004
46 * @time 10:39am
47 * @version 1.0
48 */
49
50 #include "brains/ForceManager.hpp"
51 #include "primitives/Molecule.hpp"
52 #define __OPENMD_C
53 #include "utils/simError.h"
54 #include "primitives/Bond.hpp"
55 #include "primitives/Bend.hpp"
56 #include "primitives/Torsion.hpp"
57 #include "primitives/Inversion.hpp"
58 #include "nonbonded/NonBondedInteraction.hpp"
59 #include "parallel/ForceMatrixDecomposition.hpp"
60
61 #include <cstdio>
62 #include <iostream>
63 #include <iomanip>
64
65 #include <omp.h>
66
67 using namespace std;
68 namespace OpenMD {
69
70 ForceManager::ForceManager(SimInfo * info) :
71 info_(info) {
72 forceField_ = info_->getForceField();
73 interactionMan_ = new InteractionManager();
74 fDecomp_ = new ForceMatrixDecomposition(info_, interactionMan_);
75 }
76
77 /**
78 * setupCutoffs
79 *
80 * Sets the values of cutoffRadius, switchingRadius, cutoffMethod,
81 * and cutoffPolicy
82 *
83 * cutoffRadius : realType
84 * If the cutoffRadius was explicitly set, use that value.
85 * If the cutoffRadius was not explicitly set:
86 * Are there electrostatic atoms? Use 12.0 Angstroms.
87 * No electrostatic atoms? Poll the atom types present in the
88 * simulation for suggested cutoff values (e.g. 2.5 * sigma).
89 * Use the maximum suggested value that was found.
90 *
91 * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE,
92 * or SHIFTED_POTENTIAL)
93 * If cutoffMethod was explicitly set, use that choice.
94 * If cutoffMethod was not explicitly set, use SHIFTED_FORCE
95 *
96 * cutoffPolicy : (one of MIX, MAX, TRADITIONAL)
97 * If cutoffPolicy was explicitly set, use that choice.
98 * If cutoffPolicy was not explicitly set, use TRADITIONAL
99 *
100 * switchingRadius : realType
101 * If the cutoffMethod was set to SWITCHED:
102 * If the switchingRadius was explicitly set, use that value
103 * (but do a sanity check first).
104 * If the switchingRadius was not explicitly set: use 0.85 *
105 * cutoffRadius_
106 * If the cutoffMethod was not set to SWITCHED:
107 * Set switchingRadius equal to cutoffRadius for safety.
108 */
109 void ForceManager::setupCutoffs() {
110
111 Globals* simParams_ = info_->getSimParams();
112 ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
113
114 if (simParams_->haveCutoffRadius())
115 {
116 rCut_ = simParams_->getCutoffRadius();
117 } else
118 {
119 if (info_->usesElectrostaticAtoms())
120 {
121 sprintf(painCave.errMsg, "ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n"
122 "\tOpenMD will use a default value of 12.0 angstroms"
123 "\tfor the cutoffRadius.\n");
124 painCave.isFatal = 0;
125 painCave.severity = OPENMD_INFO;
126 simError();
127 rCut_ = 12.0;
128 } else
129 {
130 RealType thisCut;
131 set<AtomType*>::iterator i;
132 set<AtomType*> atomTypes;
133 atomTypes = info_->getSimulatedAtomTypes();
134 for (i = atomTypes.begin(); i != atomTypes.end(); ++i)
135 {
136 thisCut = interactionMan_->getSuggestedCutoffRadius((*i));
137 rCut_ = max(thisCut, rCut_);
138 }
139 sprintf(painCave.errMsg, "ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n"
140 "\tOpenMD will use %lf angstroms.\n", rCut_);
141 painCave.isFatal = 0;
142 painCave.severity = OPENMD_INFO;
143 simError();
144 }
145 }
146
147 fDecomp_->setUserCutoff(rCut_);
148 interactionMan_->setCutoffRadius(rCut_);
149
150 map<string, CutoffMethod> stringToCutoffMethod;
151 stringToCutoffMethod["HARD"] = HARD;
152 stringToCutoffMethod["SWITCHED"] = SWITCHED;
153 stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;
154 stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE;
155
156 if (simParams_->haveCutoffMethod())
157 {
158 string cutMeth = toUpperCopy(simParams_->getCutoffMethod());
159 map<string, CutoffMethod>::iterator i;
160 i = stringToCutoffMethod.find(cutMeth);
161 if (i == stringToCutoffMethod.end())
162 {
163 sprintf(painCave.errMsg, "ForceManager::setupCutoffs: Could not find chosen cutoffMethod %s\n"
164 "\tShould be one of: "
165 "HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n", cutMeth.c_str());
166 painCave.isFatal = 1;
167 painCave.severity = OPENMD_ERROR;
168 simError();
169 } else
170 {
171 cutoffMethod_ = i->second;
172 }
173 } else
174 {
175 sprintf(painCave.errMsg, "ForceManager::setupCutoffs: No value was set for the cutoffMethod.\n"
176 "\tOpenMD will use SHIFTED_FORCE.\n");
177 painCave.isFatal = 0;
178 painCave.severity = OPENMD_INFO;
179 simError();
180 cutoffMethod_ = SHIFTED_FORCE;
181 }
182
183 map<string, CutoffPolicy> stringToCutoffPolicy;
184 stringToCutoffPolicy["MIX"] = MIX;
185 stringToCutoffPolicy["MAX"] = MAX;
186 stringToCutoffPolicy["TRADITIONAL"] = TRADITIONAL;
187
188 std::string cutPolicy;
189 if (forceFieldOptions_.haveCutoffPolicy())
190 {
191 cutPolicy = forceFieldOptions_.getCutoffPolicy();
192 } else if (simParams_->haveCutoffPolicy())
193 {
194 cutPolicy = simParams_->getCutoffPolicy();
195 }
196
197 if (!cutPolicy.empty())
198 {
199 toUpper(cutPolicy);
200 map<string, CutoffPolicy>::iterator i;
201 i = stringToCutoffPolicy.find(cutPolicy);
202
203 if (i == stringToCutoffPolicy.end())
204 {
205 sprintf(painCave.errMsg, "ForceManager::setupCutoffs: Could not find chosen cutoffPolicy %s\n"
206 "\tShould be one of: "
207 "MIX, MAX, or TRADITIONAL\n", cutPolicy.c_str());
208 painCave.isFatal = 1;
209 painCave.severity = OPENMD_ERROR;
210 simError();
211 } else
212 {
213 cutoffPolicy_ = i->second;
214 }
215 } else
216 {
217 sprintf(painCave.errMsg, "ForceManager::setupCutoffs: No value was set for the cutoffPolicy.\n"
218 "\tOpenMD will use TRADITIONAL.\n");
219 painCave.isFatal = 0;
220 painCave.severity = OPENMD_INFO;
221 simError();
222 cutoffPolicy_ = TRADITIONAL;
223 }
224
225 fDecomp_->setCutoffPolicy(cutoffPolicy_);
226
227 // create the switching function object:
228
229 switcher_ = new SwitchingFunction();
230
231 if (cutoffMethod_ == SWITCHED)
232 {
233 if (simParams_->haveSwitchingRadius())
234 {
235 rSwitch_ = simParams_->getSwitchingRadius();
236 if (rSwitch_ > rCut_)
237 {
238 sprintf(painCave.errMsg, "ForceManager::setupCutoffs: switchingRadius (%f) is larger "
239 "than the cutoffRadius(%f)\n", rSwitch_, rCut_);
240 painCave.isFatal = 1;
241 painCave.severity = OPENMD_ERROR;
242 simError();
243 }
244 } else
245 {
246 rSwitch_ = 0.85 * rCut_;
247 sprintf(painCave.errMsg, "ForceManager::setupCutoffs: No value was set for the switchingRadius.\n"
248 "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
249 "\tswitchingRadius = %f. for this simulation\n", rSwitch_);
250 painCave.isFatal = 0;
251 painCave.severity = OPENMD_WARNING;
252 simError();
253 }
254 } else
255 {
256 if (simParams_->haveSwitchingRadius())
257 {
258 map<string, CutoffMethod>::const_iterator it;
259 string theMeth;
260 for (it = stringToCutoffMethod.begin(); it != stringToCutoffMethod.end(); ++it)
261 {
262 if (it->second == cutoffMethod_)
263 {
264 theMeth = it->first;
265 break;
266 }
267 }
268 sprintf(painCave.errMsg, "ForceManager::setupCutoffs: the cutoffMethod (%s)\n"
269 "\tis not set to SWITCHED, so switchingRadius value\n"
270 "\twill be ignored for this simulation\n", theMeth.c_str());
271 painCave.isFatal = 0;
272 painCave.severity = OPENMD_WARNING;
273 simError();
274 }
275
276 rSwitch_ = rCut_;
277 }
278
279 // Default to cubic switching function.
280 sft_ = cubic;
281 if (simParams_->haveSwitchingFunctionType())
282 {
283 string funcType = simParams_->getSwitchingFunctionType();
284 toUpper(funcType);
285 if (funcType == "CUBIC")
286 {
287 sft_ = cubic;
288 } else
289 {
290 if (funcType == "FIFTH_ORDER_POLYNOMIAL")
291 {
292 sft_ = fifth_order_poly;
293 } else
294 {
295 // throw error
296 sprintf(painCave.errMsg,
297 "ForceManager::setupSwitching : Unknown switchingFunctionType. (Input file specified %s .)\n"
298 "\tswitchingFunctionType must be one of: "
299 "\"cubic\" or \"fifth_order_polynomial\".", funcType.c_str());
300 painCave.isFatal = 1;
301 painCave.severity = OPENMD_ERROR;
302 simError();
303 }
304 }
305 }
306 switcher_->setSwitchType(sft_);
307 switcher_->setSwitch(rSwitch_, rCut_);
308 interactionMan_->setSwitchingRadius(rSwitch_);
309 }
310
311 void ForceManager::initialize() {
312
313 if (!info_->isTopologyDone())
314 {
315
316 info_->update();
317 interactionMan_->setSimInfo(info_);
318 interactionMan_->initialize();
319
320 // We want to delay the cutoffs until after the interaction
321 // manager has set up the atom-atom interactions so that we can
322 // query them for suggested cutoff values
323 setupCutoffs();
324
325 info_->prepareTopology();
326 }
327
328 ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
329
330 // Force fields can set options on how to scale van der Waals and
331 // electrostatic interactions for atoms connected via bonds, bends
332 // and torsions in this case the topological distance between
333 // atoms is:
334 // 0 = topologically unconnected
335 // 1 = bonded together
336 // 2 = connected via a bend
337 // 3 = connected via a torsion
338
339 vdwScale_.reserve(4);
340 fill(vdwScale_.begin(), vdwScale_.end(), 0.0);
341
342 electrostaticScale_.reserve(4);
343 fill(electrostaticScale_.begin(), electrostaticScale_.end(), 0.0);
344
345 vdwScale_[0] = 1.0;
346 vdwScale_[1] = fopts.getvdw12scale();
347 vdwScale_[2] = fopts.getvdw13scale();
348 vdwScale_[3] = fopts.getvdw14scale();
349
350 electrostaticScale_[0] = 1.0;
351 electrostaticScale_[1] = fopts.getelectrostatic12scale();
352 electrostaticScale_[2] = fopts.getelectrostatic13scale();
353 electrostaticScale_[3] = fopts.getelectrostatic14scale();
354
355 fDecomp_->distributeInitialData();
356
357 initialized_ = true;
358
359 }
360
361 void ForceManager::calcForces() {
362
363 if (!initialized_)
364 initialize();
365
366 preCalculation();
367 shortRangeInteractions();
368 // longRangeInteractions();
369 // longRangeInteractionsRapaport();
370 longRangeInteractionsParallel();
371 postCalculation();
372 }
373
374 void ForceManager::preCalculation() {
375 SimInfo::MoleculeIterator mi;
376 Molecule* mol;
377 Molecule::AtomIterator ai;
378 Atom* atom;
379 Molecule::RigidBodyIterator rbIter;
380 RigidBody* rb;
381 Molecule::CutoffGroupIterator ci;
382 CutoffGroup* cg;
383
384 // forces are zeroed here, before any are accumulated.
385
386 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
387 {
388 for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai))
389 {
390 atom->zeroForcesAndTorques();
391 }
392
393 //change the positions of atoms which belong to the rigidbodies
394 for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter))
395 {
396 rb->zeroForcesAndTorques();
397 }
398
399 if (info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms())
400 {
401 for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
402 {
403 //calculate the center of mass of cutoff group
404 cg->updateCOM();
405 }
406 }
407 }
408
409 // Zero out the stress tensor
410 tau *= 0.0;
411
412 }
413
414 void ForceManager::shortRangeInteractions() {
415 Molecule* mol;
416 RigidBody* rb;
417 Bond* bond;
418 Bend* bend;
419 Torsion* torsion;
420 Inversion* inversion;
421 SimInfo::MoleculeIterator mi;
422 Molecule::RigidBodyIterator rbIter;
423 Molecule::BondIterator bondIter;
424 ;
425 Molecule::BendIterator bendIter;
426 Molecule::TorsionIterator torsionIter;
427 Molecule::InversionIterator inversionIter;
428 RealType bondPotential = 0.0;
429 RealType bendPotential = 0.0;
430 RealType torsionPotential = 0.0;
431 RealType inversionPotential = 0.0;
432
433 //calculate short range interactions
434 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
435 {
436
437 //change the positions of atoms which belong to the rigidbodies
438 for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter))
439 {
440 rb->updateAtoms();
441 }
442
443 for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter))
444 {
445 bond->calcForce();
446 bondPotential += bond->getPotential();
447 }
448
449 for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter))
450 {
451
452 RealType angle;
453 bend->calcForce(angle);
454 RealType currBendPot = bend->getPotential();
455
456 bendPotential += bend->getPotential();
457 map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
458 if (i == bendDataSets.end())
459 {
460 BendDataSet dataSet;
461 dataSet.prev.angle = dataSet.curr.angle = angle;
462 dataSet.prev.potential = dataSet.curr.potential = currBendPot;
463 dataSet.deltaV = 0.0;
464 bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend, dataSet));
465 } else
466 {
467 i->second.prev.angle = i->second.curr.angle;
468 i->second.prev.potential = i->second.curr.potential;
469 i->second.curr.angle = angle;
470 i->second.curr.potential = currBendPot;
471 i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential);
472 }
473 }
474
475 for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter))
476 {
477 RealType angle;
478 torsion->calcForce(angle);
479 RealType currTorsionPot = torsion->getPotential();
480 torsionPotential += torsion->getPotential();
481 map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
482 if (i == torsionDataSets.end())
483 {
484 TorsionDataSet dataSet;
485 dataSet.prev.angle = dataSet.curr.angle = angle;
486 dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
487 dataSet.deltaV = 0.0;
488 torsionDataSets.insert(map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
489 } else
490 {
491 i->second.prev.angle = i->second.curr.angle;
492 i->second.prev.potential = i->second.curr.potential;
493 i->second.curr.angle = angle;
494 i->second.curr.potential = currTorsionPot;
495 i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential);
496 }
497 }
498
499 for (inversion = mol->beginInversion(inversionIter); inversion != NULL; inversion = mol->nextInversion(
500 inversionIter))
501 {
502 RealType angle;
503 inversion->calcForce(angle);
504 RealType currInversionPot = inversion->getPotential();
505 inversionPotential += inversion->getPotential();
506 map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
507 if (i == inversionDataSets.end())
508 {
509 InversionDataSet dataSet;
510 dataSet.prev.angle = dataSet.curr.angle = angle;
511 dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
512 dataSet.deltaV = 0.0;
513 inversionDataSets.insert(map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
514 } else
515 {
516 i->second.prev.angle = i->second.curr.angle;
517 i->second.prev.potential = i->second.curr.potential;
518 i->second.curr.angle = angle;
519 i->second.curr.potential = currInversionPot;
520 i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential);
521 }
522 }
523 }
524
525 RealType shortRangePotential = bondPotential + bendPotential + torsionPotential + inversionPotential;
526 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
527 curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
528 curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
529 curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
530 curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
531 curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
532 }
533
534 void ForceManager::longRangeInteractionsParallel() {
535
536 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
537 DataStorage* config = &(curSnapshot->atomData);
538 DataStorage* cgConfig = &(curSnapshot->cgData);
539
540 //calculate the center of mass of cutoff group
541
542 SimInfo::MoleculeIterator mi;
543 Molecule* mol;
544 Molecule::CutoffGroupIterator ci;
545 CutoffGroup* cg;
546
547 if (info_->getNCutoffGroups() > 0)
548 {
549 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
550 {
551 for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
552 {
553 // cerr << "branch1\n";
554 // cerr << "globind = " << cg->getGlobalIndex() << ":" << __LINE__ << "\n";
555 cg->updateCOM();
556
557 // cerr << "gbI: " << cg->getGlobalIndex() << " locI: " << cg->getLocalIndex() << " x: "
558 // << cgConfig->position[cg->getLocalIndex()].x() << " y: " << cgConfig->position[cg->getLocalIndex()].y()
559 // << " z: " << cgConfig->position[cg->getLocalIndex()].z() << "\n";
560 }
561 }
562 } else
563 {
564 // center of mass of the group is the same as position of the atom
565 // if cutoff group does not exist
566 // cerr << ":" << __LINE__ << "branch2\n";
567 cgConfig->position = config->position;
568 }
569
570 fDecomp_->zeroWorkArrays();
571 fDecomp_->distributeData();
572
573 int atom1, atom2, topoDist;
574 Vector3d d_grp, dag, d;
575 RealType rgrpsq, rgrp, r2, r;
576 RealType electroMult, vdwMult;
577 RealType vij;
578 Vector3d fij, fg, f1;
579 tuple3<RealType, RealType, RealType> cuts;
580 RealType rCutSq;
581 bool in_switching_region;
582 RealType sw, dswdr, swderiv;
583 vector<int> atomListColumn, atomListRow, atomListLocal;
584
585 InteractionDataPrv idatPrv;
586
587 SelfData sdat;
588 RealType mf;
589 RealType lrPot;
590 RealType vpair;
591 potVec longRangePotential(0.0);
592 potVec workPot(0.0);
593
594 int loopStart, loopEnd;
595 sdat.pot = fDecomp_->getEmbeddingPotential();
596
597 vector<CutoffGroup *> cgs;
598
599 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
600 {
601 for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
602 {
603 cgs.push_back(cg);
604 }
605 }
606
607 loopEnd = PAIR_LOOP;
608 if (info_->requiresPrepair())
609 {
610 loopStart = PREPAIR_LOOP;
611 } else
612 {
613 loopStart = PAIR_LOOP;
614 }
615
616 for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++)
617 {
618
619 if (iLoop == loopStart)
620 {
621 bool update_nlist = fDecomp_->checkNeighborList();
622 if (update_nlist)
623 neighborMatW = fDecomp_->buildLayerBasedNeighborList();
624 }
625
626 vector<CutoffGroup *>::iterator cg1;
627 vector<CutoffGroup *>::iterator cg2;
628
629 // int nThreads = 2;
630 int chunkSize = cgs.size() / (omp_get_num_threads() * 20);
631
632 // printf("before omp loop\n");
633 #pragma omp parallel /*num_threads(nThreads)*/ default(none) shared(curSnapshot, iLoop, cgs, chunkSize) \
634 private(cg1, cg2, cuts, d_grp, rgrpsq, rCutSq, idatPrv, vij, fij, in_switching_region, dswdr, rgrp, \
635 atomListRow, atomListColumn, atom1, atom2, topoDist, d, r2, swderiv, fg, mf, dag)
636 {
637 idatPrv.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
638 idatPrv.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
639
640 // printf("Thread %d\n", omp_get_thread_num());
641 #pragma omp for schedule(dynamic, chunkSize)
642 for (cg1 = cgs.begin(); cg1 < cgs.end(); ++cg1)
643 {
644 for (cg2 = neighborMatW[(*cg1)->getGlobalIndex()].begin(); cg2 < neighborMatW[(*cg1)->getGlobalIndex()].end(); ++cg2)
645 {
646
647 cuts = fDecomp_->getGroupCutoffs((*cg1)->getGlobalIndex(), (*cg2)->getGlobalIndex());
648
649 d_grp = fDecomp_->getIntergroupVector((*cg1), (*cg2));
650 curSnapshot->wrapVector(d_grp);
651 rgrpsq = d_grp.lengthSquare();
652
653 rCutSq = cuts.second;
654
655 // printf("Thread %d\tcg1:%d\tcg2:%d d_grp\tx:%f\ty:%f\tz:%f\trgrpsq:%f\n", omp_get_thread_num(), (*cg1)->getGlobalIndex(), (*cg2)->getGlobalIndex(), d_grp.x(), d_grp.y(), d_grp.z(), rgrpsq);
656
657 if (rgrpsq < rCutSq)
658 {
659 idatPrv.rcut = cuts.first;
660 if (iLoop == PAIR_LOOP)
661 {
662 vij = 0.0;
663 fij = V3Zero;
664 }
665
666 in_switching_region = switcher_->getSwitch(rgrpsq, /*sw*/idatPrv.sw, dswdr, rgrp);
667
668 // printf("in_switching_region:%d\trgrpsq:%f\t*idatPrv.sw:%f\tdswdr:%f\trgrp:%f\n", (in_switching_region == false ? 0 : 1), rgrpsq, idatPrv.sw, dswdr, rgrp);
669
670 atomListRow = fDecomp_->getAtomsInGroupRow((*cg1)->getGlobalIndex());
671 atomListColumn = fDecomp_->getAtomsInGroupColumn((*cg2)->getGlobalIndex());
672
673 for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
674 {
675 atom1 = (*ia);
676
677 for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
678 {
679 atom2 = (*jb);
680
681 // printf("atom1:%d atom2:%d\n", atom1, atom2);
682 if (!fDecomp_->skipAtomPair(atom1, atom2))
683 {
684 idatPrv.vpair = 0.0;
685 idatPrv.pot = 0.0;
686 idatPrv.f1 = V3Zero;
687
688 fDecomp_->fillInteractionDataOMP(idatPrv, atom1, atom2);
689
690 topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
691 idatPrv.vdwMult = vdwScale_[topoDist];
692 idatPrv.electroMult = electrostaticScale_[topoDist];
693
694 // printf("topoDist:%d\tidatPrv.vdwMult:%f\tidatPrv.electroMult:%f\n", topoDist, idatPrv.vdwMult, idatPrv.electroMult);
695
696 if (atomListRow.size() == 1 && atomListColumn.size() == 1)
697 {
698 idatPrv.d = d_grp;
699 idatPrv.r2 = rgrpsq;
700 // cerr << "dgrp = " << d_grp << ":" << __LINE__ << "\n";
701 } else
702 {
703 d = fDecomp_->getInteratomicVector(atom1, atom2);
704 curSnapshot->wrapVector(d);
705 r2 = d.lengthSquare();
706 // cerr << "datm = " << d << ":" << __LINE__ << "\n";
707 idatPrv.d = d;
708 idatPrv.r2 = r2;
709 }
710
711 // printf("idatPrv.d x:%f\ty:%f\tz:%f\tidatPrv.r2:%f\n", (idatPrv.d).x(), (idatPrv.d).y(), (idatPrv.d).z(), idatPrv.r2);
712
713 // cerr << "idat.d = " << *(idat.d) << ":" << __LINE__ << "\n";
714 idatPrv.rij = sqrt((idatPrv.r2));
715 // cerr << "idat.rij = " << *(idat.rij) << "\n";
716
717 #pragma omp critical
718 {
719 interactionMan_->initializeOMP();
720 }
721
722 if (iLoop == PREPAIR_LOOP)
723 {
724 interactionMan_->doPrePairOMP(idatPrv);
725 } else
726 {
727 interactionMan_->doPairOMP(idatPrv);
728 fDecomp_->unpackInteractionDataOMP(idatPrv, atom1, atom2);
729
730 // cerr << "d = " << *(idat.d) << "\tv=" << vpair << "\tf=" << f1 << ":" << __LINE__ << "\n";
731 // printf("d x:%f y:%f z:%f vpair:%f f1 x:%f y:%f z:%f\n", idatPrv.d.x(), idatPrv.d.y(), idatPrv.d.z(), idatPrv.vpair, idatPrv.f1.x(), idatPrv.f1.y(), idatPrv.f1.z());
732 #pragma omp critical
733 {
734 vij += idatPrv.vpair;
735 fij += idatPrv.f1;
736 tau -= outProduct(idatPrv.d, idatPrv.f1);
737
738 // printf("vij:%f fij x:%f y:%f z:%f\n", vij, fij.x(), fij.y(), fij.z());
739 }
740 }
741 }
742 }
743 }
744
745 if (iLoop == PAIR_LOOP)
746 {
747 if (in_switching_region)
748 {
749 swderiv = vij * dswdr / rgrp;
750 fg = swderiv * d_grp;
751 fij += fg;
752
753 if (atomListRow.size() == 1 && atomListColumn.size() == 1)
754 {
755 #pragma omp critical
756 {
757 tau -= outProduct(idatPrv.d, fg);
758 }
759 }
760
761 for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
762 {
763 atom1 = (*ia);
764 mf = fDecomp_->getMassFactorRow(atom1);
765 // fg is the force on atom ia due to cutoff group's
766 // presence in switching region
767 fg = swderiv * d_grp * mf;
768 fDecomp_->addForceToAtomRowOMP(atom1, fg);
769
770 if (atomListRow.size() > 1)
771 {
772 if (info_->usesAtomicVirial())
773 {
774 // find the distance between the atom
775 // and the center of the cutoff group:
776 dag = fDecomp_->getAtomToGroupVectorRow(atom1, (*cg1)->getGlobalIndex());
777 #pragma omp critical
778 {
779 tau -= outProduct(dag, fg);
780 }
781 }
782 }
783 }
784 for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
785 {
786 atom2 = (*jb);
787 mf = fDecomp_->getMassFactorColumn(atom2);
788 // fg is the force on atom jb due to cutoff group's
789 // presence in switching region
790 fg = -swderiv * d_grp * mf;
791 fDecomp_->addForceToAtomColumn(atom2, fg);
792
793 if (atomListColumn.size() > 1)
794 {
795 if (info_->usesAtomicVirial())
796 {
797 // find the distance between the atom
798 // and the center of the cutoff group:
799 dag = fDecomp_->getAtomToGroupVectorColumn(atom2, (*cg2)->getGlobalIndex());
800 #pragma omp critical
801 {
802 tau -= outProduct(dag, fg);
803 }
804 }
805 }
806 }
807 }
808 //if (!SIM_uses_AtomicVirial) {
809 // tau -= outProduct(d_grp, fij);
810 //}
811 }
812 }
813 }
814 }// END: omp for loop
815 // printf("after omp loop\n");
816 }
817
818 if (iLoop == PREPAIR_LOOP)
819 {
820 if (info_->requiresPrepair())
821 {
822
823 fDecomp_->collectIntermediateData();
824
825 for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
826 {
827 fDecomp_->fillSelfData(sdat, atom1);
828 interactionMan_->doPreForce(sdat);
829 }
830
831 fDecomp_->distributeIntermediateData();
832
833 }
834 }
835 }
836
837 fDecomp_->collectData();
838
839 if (info_->requiresSelfCorrection())
840 {
841
842 for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
843 {
844 fDecomp_->fillSelfData(sdat, atom1);
845 interactionMan_->doSelfCorrection(sdat);
846 }
847
848 }
849
850 longRangePotential = *(fDecomp_->getEmbeddingPotential()) + *(fDecomp_->getPairwisePotential());
851
852 lrPot = longRangePotential.sum();
853
854 //store the tau and long range potential
855 curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
856 curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
857 curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
858 }
859
860 void ForceManager::longRangeInteractionsRapaport() {
861
862 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
863 DataStorage* config = &(curSnapshot->atomData);
864 DataStorage* cgConfig = &(curSnapshot->cgData);
865
866 //calculate the center of mass of cutoff group
867
868 SimInfo::MoleculeIterator mi;
869 Molecule* mol;
870 Molecule::CutoffGroupIterator ci;
871 CutoffGroup* cg;
872
873 if (info_->getNCutoffGroups() > 0)
874 {
875 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
876 {
877 for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
878 {
879 // cerr << "branch1\n";
880 // cerr << "globind = " << cg->getGlobalIndex() << ":" << __LINE__ << "\n";
881 cg->updateCOM();
882
883 // cerr << "gbI: " << cg->getGlobalIndex() << " locI: " << cg->getLocalIndex() << " x: "
884 // << cgConfig->position[cg->getLocalIndex()].x() << " y: " << cgConfig->position[cg->getLocalIndex()].y()
885 // << " z: " << cgConfig->position[cg->getLocalIndex()].z() << "\n";
886 }
887 }
888 } else
889 {
890 // center of mass of the group is the same as position of the atom
891 // if cutoff group does not exist
892 // cerr << ":" << __LINE__ << "branch2\n";
893 cgConfig->position = config->position;
894 }
895
896 fDecomp_->zeroWorkArrays();
897 fDecomp_->distributeData();
898
899 int atom1, atom2, topoDist;
900 CutoffGroup *cg1;
901 Vector3d d_grp, dag, d;
902 RealType rgrpsq, rgrp, r2, r;
903 RealType electroMult, vdwMult;
904 RealType vij;
905 Vector3d fij, fg, f1;
906 tuple3<RealType, RealType, RealType> cuts;
907 RealType rCutSq;
908 bool in_switching_region;
909 RealType sw, dswdr, swderiv;
910 vector<int> atomListColumn, atomListRow, atomListLocal;
911 InteractionData idat;
912 SelfData sdat;
913 RealType mf;
914 RealType lrPot;
915 RealType vpair;
916 potVec longRangePotential(0.0);
917 potVec workPot(0.0);
918
919 int loopStart, loopEnd;
920
921 idat.vdwMult = &vdwMult;
922 idat.electroMult = &electroMult;
923 idat.pot = &workPot;
924 sdat.pot = fDecomp_->getEmbeddingPotential();
925 idat.vpair = &vpair;
926 idat.f1 = &f1;
927 idat.sw = &sw;
928 idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
929 idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
930
931 loopEnd = PAIR_LOOP;
932 if (info_->requiresPrepair())
933 {
934 loopStart = PREPAIR_LOOP;
935 } else
936 {
937 loopStart = PAIR_LOOP;
938 }
939
940 for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++)
941 {
942
943 if (iLoop == loopStart)
944 {
945 bool update_nlist = fDecomp_->checkNeighborList();
946 if (update_nlist)
947 neighborMatW = fDecomp_->buildLayerBasedNeighborList();
948 }
949
950 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
951 {
952 for (cg1 = mol->beginCutoffGroup(ci); cg1 != NULL; cg1 = mol->nextCutoffGroup(ci))
953 {
954 // printf("Thread %d executes loop iteration %d\n", omp_get_thread_num(), i);
955 for (vector<CutoffGroup *>::iterator cg2 = neighborMatW[cg1->getGlobalIndex()].begin(); cg2
956 != neighborMatW[cg1->getGlobalIndex()].end(); ++cg2)
957 {
958
959 cuts = fDecomp_->getGroupCutoffs(cg1->getGlobalIndex(), (*cg2)->getGlobalIndex());
960
961 d_grp = fDecomp_->getIntergroupVector(cg1, (*cg2));
962 curSnapshot->wrapVector(d_grp);
963 rgrpsq = d_grp.lengthSquare();
964
965 rCutSq = cuts.second;
966
967 if (rgrpsq < rCutSq)
968 {
969 idat.rcut = &cuts.first;
970 if (iLoop == PAIR_LOOP)
971 {
972 vij = 0.0;
973 fij = V3Zero;
974 }
975
976 in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, rgrp);
977
978 atomListRow = fDecomp_->getAtomsInGroupRow(cg1->getGlobalIndex());
979 atomListColumn = fDecomp_->getAtomsInGroupColumn((*cg2)->getGlobalIndex());
980
981 for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
982 {
983 atom1 = (*ia);
984
985 for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
986 {
987 atom2 = (*jb);
988
989 if (!fDecomp_->skipAtomPair(atom1, atom2))
990 {
991 vpair = 0.0;
992 workPot = 0.0;
993 f1 = V3Zero;
994
995 fDecomp_->fillInteractionData(idat, atom1, atom2);
996
997 topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
998 vdwMult = vdwScale_[topoDist];
999 electroMult = electrostaticScale_[topoDist];
1000
1001 if (atomListRow.size() == 1 && atomListColumn.size() == 1)
1002 {
1003 idat.d = &d_grp;
1004 idat.r2 = &rgrpsq;
1005 // cerr << "dgrp = " << d_grp << ":" << __LINE__ << "\n";
1006 } else
1007 {
1008 d = fDecomp_->getInteratomicVector(atom1, atom2);
1009 curSnapshot->wrapVector(d);
1010 r2 = d.lengthSquare();
1011 // cerr << "datm = " << d << ":" << __LINE__ << "\n";
1012 idat.d = &d;
1013 idat.r2 = &r2;
1014 }
1015
1016 // cerr << "idat.d = " << *(idat.d) << ":" << __LINE__ << "\n";
1017 r = sqrt(*(idat.r2));
1018 idat.rij = &r;
1019 // cerr << "idat.rij = " << *(idat.rij) << "\n";
1020
1021 if (iLoop == PREPAIR_LOOP)
1022 {
1023 interactionMan_->doPrePair(idat);
1024 } else
1025 {
1026 interactionMan_->doPair(idat);
1027 fDecomp_->unpackInteractionData(idat, atom1, atom2);
1028
1029 // cerr << "d = " << *(idat.d) << "\tv=" << vpair << "\tf=" << f1 << ":" << __LINE__ << "\n";
1030 vij += vpair;
1031 fij += f1;
1032 tau -= outProduct(*(idat.d), f1);
1033 }
1034 }
1035 }
1036 }
1037
1038 if (iLoop == PAIR_LOOP)
1039 {
1040 if (in_switching_region)
1041 {
1042 swderiv = vij * dswdr / rgrp;
1043 fg = swderiv * d_grp;
1044 fij += fg;
1045
1046 if (atomListRow.size() == 1 && atomListColumn.size() == 1)
1047 {
1048 tau -= outProduct(*(idat.d), fg);
1049 }
1050
1051 for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
1052 {
1053 atom1 = (*ia);
1054 mf = fDecomp_->getMassFactorRow(atom1);
1055 // fg is the force on atom ia due to cutoff group's
1056 // presence in switching region
1057 fg = swderiv * d_grp * mf;
1058 fDecomp_->addForceToAtomRow(atom1, fg);
1059
1060 if (atomListRow.size() > 1)
1061 {
1062 if (info_->usesAtomicVirial())
1063 {
1064 // find the distance between the atom
1065 // and the center of the cutoff group:
1066 dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1->getGlobalIndex());
1067 tau -= outProduct(dag, fg);
1068 }
1069 }
1070 }
1071 for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
1072 {
1073 atom2 = (*jb);
1074 mf = fDecomp_->getMassFactorColumn(atom2);
1075 // fg is the force on atom jb due to cutoff group's
1076 // presence in switching region
1077 fg = -swderiv * d_grp * mf;
1078 fDecomp_->addForceToAtomColumn(atom2, fg);
1079
1080 if (atomListColumn.size() > 1)
1081 {
1082 if (info_->usesAtomicVirial())
1083 {
1084 // find the distance between the atom
1085 // and the center of the cutoff group:
1086 dag = fDecomp_->getAtomToGroupVectorColumn(atom2, (*cg2)->getGlobalIndex());
1087 tau -= outProduct(dag, fg);
1088 }
1089 }
1090 }
1091 }
1092 //if (!SIM_uses_AtomicVirial) {
1093 // tau -= outProduct(d_grp, fij);
1094 //}
1095 }
1096 }
1097 }
1098 }
1099 }
1100
1101 if (iLoop == PREPAIR_LOOP)
1102 {
1103 if (info_->requiresPrepair())
1104 {
1105
1106 fDecomp_->collectIntermediateData();
1107
1108 for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
1109 {
1110 fDecomp_->fillSelfData(sdat, atom1);
1111 interactionMan_->doPreForce(sdat);
1112 }
1113
1114 fDecomp_->distributeIntermediateData();
1115
1116 }
1117 }
1118 }
1119
1120 fDecomp_->collectData();
1121
1122 if (info_->requiresSelfCorrection())
1123 {
1124
1125 for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
1126 {
1127 fDecomp_->fillSelfData(sdat, atom1);
1128 interactionMan_->doSelfCorrection(sdat);
1129 }
1130
1131 }
1132
1133 longRangePotential = *(fDecomp_->getEmbeddingPotential()) + *(fDecomp_->getPairwisePotential());
1134
1135 lrPot = longRangePotential.sum();
1136
1137 //store the tau and long range potential
1138 curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
1139 curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
1140 curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
1141 }
1142
1143 void ForceManager::longRangeInteractions() {
1144
1145 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
1146 DataStorage* config = &(curSnapshot->atomData);
1147 DataStorage* cgConfig = &(curSnapshot->cgData);
1148
1149 //calculate the center of mass of cutoff group
1150
1151 SimInfo::MoleculeIterator mi;
1152 Molecule* mol;
1153 Molecule::CutoffGroupIterator ci;
1154 CutoffGroup* cg;
1155
1156 if (info_->getNCutoffGroups() > 0)
1157 {
1158 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
1159 {
1160 for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
1161 {
1162 cerr << "branch1\n";
1163 cerr << "globind = " << cg->getGlobalIndex() << "\n";
1164 cg->updateCOM();
1165 }
1166 }
1167 } else
1168 {
1169 // center of mass of the group is the same as position of the atom
1170 // if cutoff group does not exist
1171 cerr << "branch2\n";
1172 cgConfig->position = config->position;
1173 }
1174
1175 fDecomp_->zeroWorkArrays();
1176 fDecomp_->distributeData();
1177
1178 int cg1, cg2, atom1, atom2, topoDist;
1179 Vector3d d_grp, dag, d;
1180 RealType rgrpsq, rgrp, r2, r;
1181 RealType electroMult, vdwMult;
1182 RealType vij;
1183 Vector3d fij, fg, f1;
1184 tuple3<RealType, RealType, RealType> cuts;
1185 RealType rCutSq;
1186 bool in_switching_region;
1187 RealType sw, dswdr, swderiv;
1188 vector<int> atomListColumn, atomListRow, atomListLocal;
1189 InteractionData idat;
1190 SelfData sdat;
1191 RealType mf;
1192 RealType lrPot;
1193 RealType vpair;
1194 potVec longRangePotential(0.0);
1195 potVec workPot(0.0);
1196
1197 int loopStart, loopEnd;
1198
1199 idat.vdwMult = &vdwMult;
1200 idat.electroMult = &electroMult;
1201 idat.pot = &workPot;
1202 sdat.pot = fDecomp_->getEmbeddingPotential();
1203 idat.vpair = &vpair;
1204 idat.f1 = &f1;
1205 idat.sw = &sw;
1206 idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
1207 idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
1208
1209 loopEnd = PAIR_LOOP;
1210 if (info_->requiresPrepair())
1211 {
1212 loopStart = PREPAIR_LOOP;
1213 } else
1214 {
1215 loopStart = PAIR_LOOP;
1216 }
1217
1218 for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++)
1219 {
1220
1221 if (iLoop == loopStart)
1222 {
1223 bool update_nlist = fDecomp_->checkNeighborList();
1224 if (update_nlist)
1225 neighborList = fDecomp_->buildNeighborList();
1226
1227 }
1228
1229 for (vector<pair<int, int> >::iterator it = neighborList.begin(); it != neighborList.end(); ++it)
1230 {
1231 cg1 = (*it).first;
1232 cg2 = (*it).second;
1233
1234 cuts = fDecomp_->getGroupCutoffs(cg1, cg2);
1235
1236 d_grp = fDecomp_->getIntergroupVector(cg1, cg2);
1237 curSnapshot->wrapVector(d_grp);
1238 rgrpsq = d_grp.lengthSquare();
1239
1240 rCutSq = cuts.second;
1241
1242 if (rgrpsq < rCutSq)
1243 {
1244 idat.rcut = &cuts.first;
1245 if (iLoop == PAIR_LOOP)
1246 {
1247 vij = 0.0;
1248 fij = V3Zero;
1249 }
1250
1251 in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, rgrp);
1252
1253 atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
1254 atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
1255
1256 for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
1257 {
1258 atom1 = (*ia);
1259
1260 for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
1261 {
1262 atom2 = (*jb);
1263
1264 if (!fDecomp_->skipAtomPair(atom1, atom2))
1265 {
1266 vpair = 0.0;
1267 workPot = 0.0;
1268 f1 = V3Zero;
1269
1270 fDecomp_->fillInteractionData(idat, atom1, atom2);
1271
1272 topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
1273 vdwMult = vdwScale_[topoDist];
1274 electroMult = electrostaticScale_[topoDist];
1275
1276 if (atomListRow.size() == 1 && atomListColumn.size() == 1)
1277 {
1278 idat.d = &d_grp;
1279 idat.r2 = &rgrpsq;
1280 cerr << "dgrp = " << d_grp << "\n";
1281 } else
1282 {
1283 d = fDecomp_->getInteratomicVector(atom1, atom2);
1284 curSnapshot->wrapVector(d);
1285 r2 = d.lengthSquare();
1286 cerr << "datm = " << d << "\n";
1287 idat.d = &d;
1288 idat.r2 = &r2;
1289 }
1290
1291 cerr << "idat.d = " << *(idat.d) << "\n";
1292 r = sqrt(*(idat.r2));
1293 idat.rij = &r;
1294
1295 if (iLoop == PREPAIR_LOOP)
1296 {
1297 interactionMan_->doPrePair(idat);
1298 } else
1299 {
1300 interactionMan_->doPair(idat);
1301 fDecomp_->unpackInteractionData(idat, atom1, atom2);
1302
1303 cerr << "d = " << *(idat.d) << "\tv=" << vpair << "\tf=" << f1 << "\n";
1304 vij += vpair;
1305 fij += f1;
1306 tau -= outProduct(*(idat.d), f1);
1307 }
1308 }
1309 }
1310 }
1311
1312 if (iLoop == PAIR_LOOP)
1313 {
1314 if (in_switching_region)
1315 {
1316 swderiv = vij * dswdr / rgrp;
1317 fg = swderiv * d_grp;
1318 fij += fg;
1319
1320 if (atomListRow.size() == 1 && atomListColumn.size() == 1)
1321 {
1322 tau -= outProduct(*(idat.d), fg);
1323 }
1324
1325 for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
1326 {
1327 atom1 = (*ia);
1328 mf = fDecomp_->getMassFactorRow(atom1);
1329 // fg is the force on atom ia due to cutoff group's
1330 // presence in switching region
1331 fg = swderiv * d_grp * mf;
1332 fDecomp_->addForceToAtomRow(atom1, fg);
1333
1334 if (atomListRow.size() > 1)
1335 {
1336 if (info_->usesAtomicVirial())
1337 {
1338 // find the distance between the atom
1339 // and the center of the cutoff group:
1340 dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
1341 tau -= outProduct(dag, fg);
1342 }
1343 }
1344 }
1345 for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
1346 {
1347 atom2 = (*jb);
1348 mf = fDecomp_->getMassFactorColumn(atom2);
1349 // fg is the force on atom jb due to cutoff group's
1350 // presence in switching region
1351 fg = -swderiv * d_grp * mf;
1352 fDecomp_->addForceToAtomColumn(atom2, fg);
1353
1354 if (atomListColumn.size() > 1)
1355 {
1356 if (info_->usesAtomicVirial())
1357 {
1358 // find the distance between the atom
1359 // and the center of the cutoff group:
1360 dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
1361 tau -= outProduct(dag, fg);
1362 }
1363 }
1364 }
1365 }
1366 //if (!SIM_uses_AtomicVirial) {
1367 // tau -= outProduct(d_grp, fij);
1368 //}
1369 }
1370 }
1371 }
1372
1373 if (iLoop == PREPAIR_LOOP)
1374 {
1375 if (info_->requiresPrepair())
1376 {
1377
1378 fDecomp_->collectIntermediateData();
1379
1380 for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
1381 {
1382 fDecomp_->fillSelfData(sdat, atom1);
1383 interactionMan_->doPreForce(sdat);
1384 }
1385
1386 fDecomp_->distributeIntermediateData();
1387
1388 }
1389 }
1390
1391 }
1392
1393 fDecomp_->collectData();
1394
1395 if (info_->requiresSelfCorrection())
1396 {
1397
1398 for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
1399 {
1400 fDecomp_->fillSelfData(sdat, atom1);
1401 interactionMan_->doSelfCorrection(sdat);
1402 }
1403
1404 }
1405
1406 longRangePotential = *(fDecomp_->getEmbeddingPotential()) + *(fDecomp_->getPairwisePotential());
1407
1408 lrPot = longRangePotential.sum();
1409
1410 //store the tau and long range potential
1411 curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
1412 curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
1413 curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
1414 }
1415
1416 void ForceManager::postCalculation() {
1417 SimInfo::MoleculeIterator mi;
1418 Molecule* mol;
1419 Molecule::RigidBodyIterator rbIter;
1420 RigidBody* rb;
1421 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
1422
1423 // collect the atomic forces onto rigid bodies
1424
1425 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
1426 {
1427 for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter))
1428 {
1429 Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
1430 tau += rbTau;
1431 }
1432 }
1433
1434 #ifdef IS_MPI
1435 Mat3x3d tmpTau(tau);
1436 MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
1437 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1438 #endif
1439 curSnapshot->statData.setTau(tau);
1440 }
1441
1442 } //end namespace OpenMD

Properties

Name Value
svn:executable *
svn:keywords Author Id Revision Date