ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/devel_omp/src/brains/ForceManager.cpp
Revision: 1598
Committed: Wed Jul 27 14:26:53 2011 UTC (14 years, 3 months ago) by mciznick
File size: 33040 byte(s)
Log Message: 
Updated ordered neighbor list generation.

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 postCalculation();
371 }
372
373 void ForceManager::preCalculation() {
374 SimInfo::MoleculeIterator mi;
375 Molecule* mol;
376 Molecule::AtomIterator ai;
377 Atom* atom;
378 Molecule::RigidBodyIterator rbIter;
379 RigidBody* rb;
380 Molecule::CutoffGroupIterator ci;
381 CutoffGroup* cg;
382
383 // forces are zeroed here, before any are accumulated.
384
385 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
386 {
387 for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai))
388 {
389 atom->zeroForcesAndTorques();
390 }
391
392 //change the positions of atoms which belong to the rigidbodies
393 for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter))
394 {
395 rb->zeroForcesAndTorques();
396 }
397
398 if (info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms())
399 {
400 for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
401 {
402 //calculate the center of mass of cutoff group
403 cg->updateCOM();
404 }
405 }
406 }
407
408 // Zero out the stress tensor
409 tau *= 0.0;
410
411 }
412
413 void ForceManager::shortRangeInteractions() {
414 Molecule* mol;
415 RigidBody* rb;
416 Bond* bond;
417 Bend* bend;
418 Torsion* torsion;
419 Inversion* inversion;
420 SimInfo::MoleculeIterator mi;
421 Molecule::RigidBodyIterator rbIter;
422 Molecule::BondIterator bondIter;
423 ;
424 Molecule::BendIterator bendIter;
425 Molecule::TorsionIterator torsionIter;
426 Molecule::InversionIterator inversionIter;
427 RealType bondPotential = 0.0;
428 RealType bendPotential = 0.0;
429 RealType torsionPotential = 0.0;
430 RealType inversionPotential = 0.0;
431
432 //calculate short range interactions
433 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
434 {
435
436 //change the positions of atoms which belong to the rigidbodies
437 for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter))
438 {
439 rb->updateAtoms();
440 }
441
442 for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter))
443 {
444 bond->calcForce();
445 bondPotential += bond->getPotential();
446 }
447
448 for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter))
449 {
450
451 RealType angle;
452 bend->calcForce(angle);
453 RealType currBendPot = bend->getPotential();
454
455 bendPotential += bend->getPotential();
456 map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
457 if (i == bendDataSets.end())
458 {
459 BendDataSet dataSet;
460 dataSet.prev.angle = dataSet.curr.angle = angle;
461 dataSet.prev.potential = dataSet.curr.potential = currBendPot;
462 dataSet.deltaV = 0.0;
463 bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend, dataSet));
464 } else
465 {
466 i->second.prev.angle = i->second.curr.angle;
467 i->second.prev.potential = i->second.curr.potential;
468 i->second.curr.angle = angle;
469 i->second.curr.potential = currBendPot;
470 i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential);
471 }
472 }
473
474 for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter))
475 {
476 RealType angle;
477 torsion->calcForce(angle);
478 RealType currTorsionPot = torsion->getPotential();
479 torsionPotential += torsion->getPotential();
480 map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
481 if (i == torsionDataSets.end())
482 {
483 TorsionDataSet dataSet;
484 dataSet.prev.angle = dataSet.curr.angle = angle;
485 dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
486 dataSet.deltaV = 0.0;
487 torsionDataSets.insert(map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
488 } else
489 {
490 i->second.prev.angle = i->second.curr.angle;
491 i->second.prev.potential = i->second.curr.potential;
492 i->second.curr.angle = angle;
493 i->second.curr.potential = currTorsionPot;
494 i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential);
495 }
496 }
497
498 for (inversion = mol->beginInversion(inversionIter); inversion != NULL; inversion = mol->nextInversion(
499 inversionIter))
500 {
501 RealType angle;
502 inversion->calcForce(angle);
503 RealType currInversionPot = inversion->getPotential();
504 inversionPotential += inversion->getPotential();
505 map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
506 if (i == inversionDataSets.end())
507 {
508 InversionDataSet dataSet;
509 dataSet.prev.angle = dataSet.curr.angle = angle;
510 dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
511 dataSet.deltaV = 0.0;
512 inversionDataSets.insert(map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
513 } else
514 {
515 i->second.prev.angle = i->second.curr.angle;
516 i->second.prev.potential = i->second.curr.potential;
517 i->second.curr.angle = angle;
518 i->second.curr.potential = currInversionPot;
519 i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential);
520 }
521 }
522 }
523
524 RealType shortRangePotential = bondPotential + bendPotential + torsionPotential + inversionPotential;
525 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
526 curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
527 curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
528 curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
529 curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
530 curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
531 }
532
533 void ForceManager::longRangeInteractionsRapaport() {
534
535 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
536 DataStorage* config = &(curSnapshot->atomData);
537 DataStorage* cgConfig = &(curSnapshot->cgData);
538
539 //calculate the center of mass of cutoff group
540
541 SimInfo::MoleculeIterator mi;
542 Molecule* mol;
543 Molecule::CutoffGroupIterator ci;
544 CutoffGroup* cg;
545
546 if (info_->getNCutoffGroups() > 0)
547 {
548 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
549 {
550 for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
551 {
552 // cerr << "branch1\n";
553 // cerr << "globind = " << cg->getGlobalIndex() << ":" << __LINE__ << "\n";
554 cg->updateCOM();
555
556 // cerr << "gbI: " << cg->getGlobalIndex() << " locI: " << cg->getLocalIndex() << " x: "
557 // << cgConfig->position[cg->getLocalIndex()].x() << " y: " << cgConfig->position[cg->getLocalIndex()].y()
558 // << " z: " << cgConfig->position[cg->getLocalIndex()].z() << "\n";
559 }
560 }
561 } else
562 {
563 // center of mass of the group is the same as position of the atom
564 // if cutoff group does not exist
565 // cerr << ":" << __LINE__ << "branch2\n";
566 cgConfig->position = config->position;
567 }
568
569 fDecomp_->zeroWorkArrays();
570 fDecomp_->distributeData();
571
572 int atom1, atom2, topoDist;
573 CutoffGroup *cg1;
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 InteractionData idat;
585 SelfData sdat;
586 RealType mf;
587 RealType lrPot;
588 RealType vpair;
589 potVec longRangePotential(0.0);
590 potVec workPot(0.0);
591
592 int loopStart, loopEnd;
593
594 idat.vdwMult = &vdwMult;
595 idat.electroMult = &electroMult;
596 idat.pot = &workPot;
597 sdat.pot = fDecomp_->getEmbeddingPotential();
598 idat.vpair = &vpair;
599 idat.f1 = &f1;
600 idat.sw = &sw;
601 idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
602 idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
603
604 loopEnd = PAIR_LOOP;
605 if (info_->requiresPrepair())
606 {
607 loopStart = PREPAIR_LOOP;
608 } else
609 {
610 loopStart = PAIR_LOOP;
611 }
612
613 for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++)
614 {
615
616 if (iLoop == loopStart)
617 {
618 bool update_nlist = fDecomp_->checkNeighborList();
619 if (update_nlist)
620 neighborMatW = fDecomp_->buildLayerBasedNeighborList();
621 }
622
623 // printf("before omp loop\n");
624 //#pragma omp parallel for num_threads(3) default(none) shared(curSnapshot, idat, iLoop, sw, cerr) \
625 private(i, j, cg1, cg2, cuts, d_grp, rgrpsq, rCutSq, vij, fij, in_switching_region, rgrp, dswdr, atomListRow, atomListColumn, atom1, atom2, vpair, workPot, f1, topoDist, vdwMult, electroMult, d, r2, r, swderiv, fg, mf, dag)
626 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
627 {
628 for (cg1 = mol->beginCutoffGroup(ci); cg1 != NULL; cg1 = mol->nextCutoffGroup(ci))
629 {
630 // printf("Thread %d executes loop iteration %d\n", omp_get_thread_num(), i);
631 for (vector<CutoffGroup *>::iterator cg2 = neighborMatW[cg1->getGlobalIndex()].begin(); cg2 != neighborMatW[cg1->getGlobalIndex()].end(); ++cg2)
632 {
633
634 cuts = fDecomp_->getGroupCutoffs(cg1->getGlobalIndex(), (*cg2)->getGlobalIndex());
635
636 d_grp = fDecomp_->getIntergroupVector(cg1, (*cg2));
637 curSnapshot->wrapVector(d_grp);
638 rgrpsq = d_grp.lengthSquare();
639
640 rCutSq = cuts.second;
641
642 if (rgrpsq < rCutSq)
643 {
644 idat.rcut = &cuts.first;
645 if (iLoop == PAIR_LOOP)
646 {
647 vij = 0.0;
648 fij = V3Zero;
649 }
650
651 in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, rgrp);
652
653 atomListRow = fDecomp_->getAtomsInGroupRow(cg1->getGlobalIndex());
654 atomListColumn = fDecomp_->getAtomsInGroupColumn((*cg2)->getGlobalIndex());
655
656 for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
657 {
658 atom1 = (*ia);
659
660 for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
661 {
662 atom2 = (*jb);
663
664 if (!fDecomp_->skipAtomPair(atom1, atom2))
665 {
666 vpair = 0.0;
667 workPot = 0.0;
668 f1 = V3Zero;
669
670 fDecomp_->fillInteractionData(idat, atom1, atom2);
671
672 topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
673 vdwMult = vdwScale_[topoDist];
674 electroMult = electrostaticScale_[topoDist];
675
676 if (atomListRow.size() == 1 && atomListColumn.size() == 1)
677 {
678 idat.d = &d_grp;
679 idat.r2 = &rgrpsq;
680 // cerr << "dgrp = " << d_grp << ":" << __LINE__ << "\n";
681 } else
682 {
683 d = fDecomp_->getInteratomicVector(atom1, atom2);
684 curSnapshot->wrapVector(d);
685 r2 = d.lengthSquare();
686 // cerr << "datm = " << d << ":" << __LINE__ << "\n";
687 idat.d = &d;
688 idat.r2 = &r2;
689 }
690
691 // cerr << "idat.d = " << *(idat.d) << ":" << __LINE__ << "\n";
692 r = sqrt(*(idat.r2));
693 idat.rij = &r;
694 // cerr << "idat.rij = " << *(idat.rij) << "\n";
695
696 if (iLoop == PREPAIR_LOOP)
697 {
698 interactionMan_->doPrePair(idat);
699 } else
700 {
701 interactionMan_->doPair(idat);
702 fDecomp_->unpackInteractionData(idat, atom1, atom2);
703
704 // cerr << "d = " << *(idat.d) << "\tv=" << vpair << "\tf=" << f1 << ":" << __LINE__ << "\n";
705 vij += vpair;
706 fij += f1;
707 tau -= outProduct(*(idat.d), f1);
708 }
709 }
710 }
711 }
712
713 if (iLoop == PAIR_LOOP)
714 {
715 if (in_switching_region)
716 {
717 swderiv = vij * dswdr / rgrp;
718 fg = swderiv * d_grp;
719 fij += fg;
720
721 if (atomListRow.size() == 1 && atomListColumn.size() == 1)
722 {
723 tau -= outProduct(*(idat.d), fg);
724 }
725
726 for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
727 {
728 atom1 = (*ia);
729 mf = fDecomp_->getMassFactorRow(atom1);
730 // fg is the force on atom ia due to cutoff group's
731 // presence in switching region
732 fg = swderiv * d_grp * mf;
733 fDecomp_->addForceToAtomRow(atom1, fg);
734
735 if (atomListRow.size() > 1)
736 {
737 if (info_->usesAtomicVirial())
738 {
739 // find the distance between the atom
740 // and the center of the cutoff group:
741 dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1->getGlobalIndex());
742 tau -= outProduct(dag, fg);
743 }
744 }
745 }
746 for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
747 {
748 atom2 = (*jb);
749 mf = fDecomp_->getMassFactorColumn(atom2);
750 // fg is the force on atom jb due to cutoff group's
751 // presence in switching region
752 fg = -swderiv * d_grp * mf;
753 fDecomp_->addForceToAtomColumn(atom2, fg);
754
755 if (atomListColumn.size() > 1)
756 {
757 if (info_->usesAtomicVirial())
758 {
759 // find the distance between the atom
760 // and the center of the cutoff group:
761 dag = fDecomp_->getAtomToGroupVectorColumn(atom2, (*cg2)->getGlobalIndex());
762 tau -= outProduct(dag, fg);
763 }
764 }
765 }
766 }
767 //if (!SIM_uses_AtomicVirial) {
768 // tau -= outProduct(d_grp, fij);
769 //}
770 }
771 }
772 }
773 }
774 }// END: omp for loop
775 // printf("after omp loop\n");
776
777 if (iLoop == PREPAIR_LOOP)
778 {
779 if (info_->requiresPrepair())
780 {
781
782 fDecomp_->collectIntermediateData();
783
784 for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
785 {
786 fDecomp_->fillSelfData(sdat, atom1);
787 interactionMan_->doPreForce(sdat);
788 }
789
790 fDecomp_->distributeIntermediateData();
791
792 }
793 }
794 }
795
796 fDecomp_->collectData();
797
798 if (info_->requiresSelfCorrection())
799 {
800
801 for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
802 {
803 fDecomp_->fillSelfData(sdat, atom1);
804 interactionMan_->doSelfCorrection(sdat);
805 }
806
807 }
808
809 longRangePotential = *(fDecomp_->getEmbeddingPotential()) + *(fDecomp_->getPairwisePotential());
810
811 lrPot = longRangePotential.sum();
812
813 //store the tau and long range potential
814 curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
815 curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
816 curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
817 }
818
819 void ForceManager::longRangeInteractions() {
820
821 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
822 DataStorage* config = &(curSnapshot->atomData);
823 DataStorage* cgConfig = &(curSnapshot->cgData);
824
825 //calculate the center of mass of cutoff group
826
827 SimInfo::MoleculeIterator mi;
828 Molecule* mol;
829 Molecule::CutoffGroupIterator ci;
830 CutoffGroup* cg;
831
832 if (info_->getNCutoffGroups() > 0)
833 {
834 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
835 {
836 for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
837 {
838 cerr << "branch1\n";
839 cerr << "globind = " << cg->getGlobalIndex() << "\n";
840 cg->updateCOM();
841 }
842 }
843 } else
844 {
845 // center of mass of the group is the same as position of the atom
846 // if cutoff group does not exist
847 cerr << "branch2\n";
848 cgConfig->position = config->position;
849 }
850
851 fDecomp_->zeroWorkArrays();
852 fDecomp_->distributeData();
853
854 int cg1, cg2, atom1, atom2, topoDist;
855 Vector3d d_grp, dag, d;
856 RealType rgrpsq, rgrp, r2, r;
857 RealType electroMult, vdwMult;
858 RealType vij;
859 Vector3d fij, fg, f1;
860 tuple3<RealType, RealType, RealType> cuts;
861 RealType rCutSq;
862 bool in_switching_region;
863 RealType sw, dswdr, swderiv;
864 vector<int> atomListColumn, atomListRow, atomListLocal;
865 InteractionData idat;
866 SelfData sdat;
867 RealType mf;
868 RealType lrPot;
869 RealType vpair;
870 potVec longRangePotential(0.0);
871 potVec workPot(0.0);
872
873 int loopStart, loopEnd;
874
875 idat.vdwMult = &vdwMult;
876 idat.electroMult = &electroMult;
877 idat.pot = &workPot;
878 sdat.pot = fDecomp_->getEmbeddingPotential();
879 idat.vpair = &vpair;
880 idat.f1 = &f1;
881 idat.sw = &sw;
882 idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
883 idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
884
885 loopEnd = PAIR_LOOP;
886 if (info_->requiresPrepair())
887 {
888 loopStart = PREPAIR_LOOP;
889 } else
890 {
891 loopStart = PAIR_LOOP;
892 }
893
894 for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++)
895 {
896
897 if (iLoop == loopStart)
898 {
899 bool update_nlist = fDecomp_->checkNeighborList();
900 if (update_nlist)
901 neighborList = fDecomp_->buildNeighborList();
902
903 }
904
905 for (vector<pair<int, int> >::iterator it = neighborList.begin(); it != neighborList.end(); ++it)
906 {
907 cg1 = (*it).first;
908 cg2 = (*it).second;
909
910 cuts = fDecomp_->getGroupCutoffs(cg1, cg2);
911
912 d_grp = fDecomp_->getIntergroupVector(cg1, cg2);
913 curSnapshot->wrapVector(d_grp);
914 rgrpsq = d_grp.lengthSquare();
915
916 rCutSq = cuts.second;
917
918 if (rgrpsq < rCutSq)
919 {
920 idat.rcut = &cuts.first;
921 if (iLoop == PAIR_LOOP)
922 {
923 vij = 0.0;
924 fij = V3Zero;
925 }
926
927 in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, rgrp);
928
929 atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
930 atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
931
932 for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
933 {
934 atom1 = (*ia);
935
936 for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
937 {
938 atom2 = (*jb);
939
940 if (!fDecomp_->skipAtomPair(atom1, atom2))
941 {
942 vpair = 0.0;
943 workPot = 0.0;
944 f1 = V3Zero;
945
946 fDecomp_->fillInteractionData(idat, atom1, atom2);
947
948 topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
949 vdwMult = vdwScale_[topoDist];
950 electroMult = electrostaticScale_[topoDist];
951
952 if (atomListRow.size() == 1 && atomListColumn.size() == 1)
953 {
954 idat.d = &d_grp;
955 idat.r2 = &rgrpsq;
956 cerr << "dgrp = " << d_grp << "\n";
957 } else
958 {
959 d = fDecomp_->getInteratomicVector(atom1, atom2);
960 curSnapshot->wrapVector(d);
961 r2 = d.lengthSquare();
962 cerr << "datm = " << d << "\n";
963 idat.d = &d;
964 idat.r2 = &r2;
965 }
966
967 cerr << "idat.d = " << *(idat.d) << "\n";
968 r = sqrt(*(idat.r2));
969 idat.rij = &r;
970
971 if (iLoop == PREPAIR_LOOP)
972 {
973 interactionMan_->doPrePair(idat);
974 } else
975 {
976 interactionMan_->doPair(idat);
977 fDecomp_->unpackInteractionData(idat, atom1, atom2);
978
979 cerr << "d = " << *(idat.d) << "\tv=" << vpair << "\tf=" << f1 << "\n";
980 vij += vpair;
981 fij += f1;
982 tau -= outProduct(*(idat.d), f1);
983 }
984 }
985 }
986 }
987
988 if (iLoop == PAIR_LOOP)
989 {
990 if (in_switching_region)
991 {
992 swderiv = vij * dswdr / rgrp;
993 fg = swderiv * d_grp;
994 fij += fg;
995
996 if (atomListRow.size() == 1 && atomListColumn.size() == 1)
997 {
998 tau -= outProduct(*(idat.d), fg);
999 }
1000
1001 for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
1002 {
1003 atom1 = (*ia);
1004 mf = fDecomp_->getMassFactorRow(atom1);
1005 // fg is the force on atom ia due to cutoff group's
1006 // presence in switching region
1007 fg = swderiv * d_grp * mf;
1008 fDecomp_->addForceToAtomRow(atom1, fg);
1009
1010 if (atomListRow.size() > 1)
1011 {
1012 if (info_->usesAtomicVirial())
1013 {
1014 // find the distance between the atom
1015 // and the center of the cutoff group:
1016 dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
1017 tau -= outProduct(dag, fg);
1018 }
1019 }
1020 }
1021 for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
1022 {
1023 atom2 = (*jb);
1024 mf = fDecomp_->getMassFactorColumn(atom2);
1025 // fg is the force on atom jb due to cutoff group's
1026 // presence in switching region
1027 fg = -swderiv * d_grp * mf;
1028 fDecomp_->addForceToAtomColumn(atom2, fg);
1029
1030 if (atomListColumn.size() > 1)
1031 {
1032 if (info_->usesAtomicVirial())
1033 {
1034 // find the distance between the atom
1035 // and the center of the cutoff group:
1036 dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
1037 tau -= outProduct(dag, fg);
1038 }
1039 }
1040 }
1041 }
1042 //if (!SIM_uses_AtomicVirial) {
1043 // tau -= outProduct(d_grp, fij);
1044 //}
1045 }
1046 }
1047 }
1048
1049 if (iLoop == PREPAIR_LOOP)
1050 {
1051 if (info_->requiresPrepair())
1052 {
1053
1054 fDecomp_->collectIntermediateData();
1055
1056 for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
1057 {
1058 fDecomp_->fillSelfData(sdat, atom1);
1059 interactionMan_->doPreForce(sdat);
1060 }
1061
1062 fDecomp_->distributeIntermediateData();
1063
1064 }
1065 }
1066
1067 }
1068
1069 fDecomp_->collectData();
1070
1071 if (info_->requiresSelfCorrection())
1072 {
1073
1074 for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
1075 {
1076 fDecomp_->fillSelfData(sdat, atom1);
1077 interactionMan_->doSelfCorrection(sdat);
1078 }
1079
1080 }
1081
1082 longRangePotential = *(fDecomp_->getEmbeddingPotential()) + *(fDecomp_->getPairwisePotential());
1083
1084 lrPot = longRangePotential.sum();
1085
1086 //store the tau and long range potential
1087 curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
1088 curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
1089 curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
1090 }
1091
1092 void ForceManager::postCalculation() {
1093 SimInfo::MoleculeIterator mi;
1094 Molecule* mol;
1095 Molecule::RigidBodyIterator rbIter;
1096 RigidBody* rb;
1097 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
1098
1099 // collect the atomic forces onto rigid bodies
1100
1101 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
1102 {
1103 for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter))
1104 {
1105 Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
1106 tau += rbTau;
1107 }
1108 }
1109
1110 #ifdef IS_MPI
1111 Mat3x3d tmpTau(tau);
1112 MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
1113 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1114 #endif
1115 curSnapshot->statData.setTau(tau);
1116 }
1117
1118 } //end namespace OpenMD

Properties

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