src/brains/ForceManager.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
 * [4]  Vardeman & Gezelter, in progress (2009).
 */

/**
 * @file ForceManager.cpp
 * @author tlin
 * @date 11/09/2004
 * @time 10:39am
 * @version 1.0
 */

#include "brains/ForceManager.hpp"
#include "primitives/Molecule.hpp"
#include "UseTheForce/doForces_interface.h"
#define __OPENMD_C
#include "UseTheForce/DarkSide/fInteractionMap.h"
#include "utils/simError.h"
#include "primitives/Bond.hpp"
#include "primitives/Bend.hpp"
#include "primitives/Torsion.hpp"
#include "primitives/Inversion.hpp"
namespace OpenMD {

  void ForceManager::calcForces() {

    if (!info_->isFortranInitialized()) {
      info_->update();
    }

    preCalculation();

    calcShortRangeInteraction();

    calcLongRangeInteraction();

    postCalculation();

  }

  void ForceManager::preCalculation() {
    SimInfo::MoleculeIterator mi;
    Molecule* mol;
    Molecule::AtomIterator ai;
    Atom* atom;
    Molecule::RigidBodyIterator rbIter;
    RigidBody* rb;

    // forces are zeroed here, before any are accumulated.
    // NOTE: do not rezero the forces in Fortran.

    for (mol = info_->beginMolecule(mi); mol != NULL;
         mol = info_->nextMolecule(mi)) {
      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
  atom->zeroForcesAndTorques();
      }

      //change the positions of atoms which belong to the rigidbodies
      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
           rb = mol->nextRigidBody(rbIter)) {
  rb->zeroForcesAndTorques();
      }

    }

    // Zero out the stress tensor
    tau *= 0.0;

  }

  void ForceManager::calcShortRangeInteraction() {
    Molecule* mol;
    RigidBody* rb;
    Bond* bond;
    Bend* bend;
    Torsion* torsion;
    Inversion* inversion;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    Molecule::BondIterator bondIter;;
    Molecule::BendIterator  bendIter;
    Molecule::TorsionIterator  torsionIter;
    Molecule::InversionIterator  inversionIter;
    RealType bondPotential = 0.0;
    RealType bendPotential = 0.0;
    RealType torsionPotential = 0.0;
    RealType inversionPotential = 0.0;

    //calculate short range interactions
    for (mol = info_->beginMolecule(mi); mol != NULL;
         mol = info_->nextMolecule(mi)) {

      //change the positions of atoms which belong to the rigidbodies
      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
           rb = mol->nextRigidBody(rbIter)) {
        rb->updateAtoms();
      }

      for (bond = mol->beginBond(bondIter); bond != NULL;
           bond = mol->nextBond(bondIter)) {
        bond->calcForce();
        bondPotential += bond->getPotential();
      }

      for (bend = mol->beginBend(bendIter); bend != NULL;
           bend = mol->nextBend(bendIter)) {

        RealType angle;
        bend->calcForce(angle);
        RealType currBendPot = bend->getPotential();

        bendPotential += bend->getPotential();
        std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
        if (i == bendDataSets.end()) {
          BendDataSet dataSet;
          dataSet.prev.angle = dataSet.curr.angle = angle;
          dataSet.prev.potential = dataSet.curr.potential = currBendPot;
          dataSet.deltaV = 0.0;
          bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
        }else {
          i->second.prev.angle = i->second.curr.angle;
          i->second.prev.potential = i->second.curr.potential;
          i->second.curr.angle = angle;
          i->second.curr.potential = currBendPot;
          i->second.deltaV =  fabs(i->second.curr.potential -
                                   i->second.prev.potential);
        }
      }

      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
           torsion = mol->nextTorsion(torsionIter)) {
        RealType angle;
        torsion->calcForce(angle);
        RealType currTorsionPot = torsion->getPotential();
        torsionPotential += torsion->getPotential();
        std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
        if (i == torsionDataSets.end()) {
          TorsionDataSet dataSet;
          dataSet.prev.angle = dataSet.curr.angle = angle;
          dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
          dataSet.deltaV = 0.0;
          torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
        }else {
          i->second.prev.angle = i->second.curr.angle;
          i->second.prev.potential = i->second.curr.potential;
          i->second.curr.angle = angle;
          i->second.curr.potential = currTorsionPot;
          i->second.deltaV =  fabs(i->second.curr.potential -
                                   i->second.prev.potential);
        }
      }

      for (inversion = mol->beginInversion(inversionIter);
     inversion != NULL;
           inversion = mol->nextInversion(inversionIter)) {
        RealType angle;
        inversion->calcForce(angle);
        RealType currInversionPot = inversion->getPotential();
        inversionPotential += inversion->getPotential();
        std::map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
        if (i == inversionDataSets.end()) {
          InversionDataSet dataSet;
          dataSet.prev.angle = dataSet.curr.angle = angle;
          dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
          dataSet.deltaV = 0.0;
          inversionDataSets.insert(std::map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
        }else {
          i->second.prev.angle = i->second.curr.angle;
          i->second.prev.potential = i->second.curr.potential;
          i->second.curr.angle = angle;
          i->second.curr.potential = currInversionPot;
          i->second.deltaV =  fabs(i->second.curr.potential -
                                   i->second.prev.potential);
        }
      }
    }

    RealType  shortRangePotential = bondPotential + bendPotential +
      torsionPotential +  inversionPotential;
    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
    curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
    curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
    curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
    curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;

  }

  void ForceManager::calcLongRangeInteraction() {
    Snapshot* curSnapshot;
    DataStorage* config;
    RealType* frc;
    RealType* pos;
    RealType* vel;
    RealType* trq;
    RealType* A;
    RealType* electroFrame;
    RealType* rc;
    RealType* vc;
    RealType* particlePot;

    //get current snapshot from SimInfo
    curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();

    //get array pointers
    config = &(curSnapshot->atomData);
    frc = config->getArrayPointer(DataStorage::dslForce);
    pos = config->getArrayPointer(DataStorage::dslPosition);
    vel = config->getArrayPointer(DataStorage::dslVelocity);
    trq = config->getArrayPointer(DataStorage::dslTorque);
    A   = config->getArrayPointer(DataStorage::dslAmat);
    electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
    particlePot = config->getArrayPointer(DataStorage::dslParticlePot);

    //calculate the center of mass of cutoff group
    SimInfo::MoleculeIterator mi;
    Molecule* mol;
    Molecule::CutoffGroupIterator ci;
    CutoffGroup* cg;
    Vector3d com;
    Vector3d comv;
    std::vector<Vector3d> rcGroup;
    std::vector<Vector3d> vcGroup;

    if(info_->getNCutoffGroups() > 0){

      for (mol = info_->beginMolecule(mi); mol != NULL;
           mol = info_->nextMolecule(mi)) {
        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
            cg = mol->nextCutoffGroup(ci)) {
              cg->getCOM(com);
              cg->getCOMV(comv);
              rcGroup.push_back(com);
              vcGroup.push_back(comv);
        }
      }// end for (mol)

      rc = rcGroup[0].getArrayPointer();
      vc = vcGroup[0].getArrayPointer();
    } else {
      // center of mass of the group is the same as position of the atom
      // if cutoff group does not exist
      rc = pos;
      vc = vel;
    }

    //initialize data before passing to fortran
    RealType longRangePotential[LR_POT_TYPES];
    RealType lrPot = 0.0;
    Vector3d totalDipole;
    Jv_ = 0.0;

    
    int isError = 0;

    for (int i=0; i<LR_POT_TYPES;i++){
      longRangePotential[i]=0.0; //Initialize array
    }

    doForceLoop(pos,
                vel,
                rc,
                vc,
                A,
                electroFrame,
                frc,
                trq,
                tau.getArrayPointer(),
                Jv_.getArrayPointer(),
                longRangePotential,
                particlePot,
                &isError );

    if( isError ){
      sprintf( painCave.errMsg,
         "Error returned from the fortran force calculation.\n" );
      painCave.isFatal = 1;
      simError();
    }
    for (int i=0; i<LR_POT_TYPES;i++){
      lrPot += longRangePotential[i]; //Quick hack
    }

    // grab the simulation box dipole moment if specified
    if (info_->getCalcBoxDipole()){
      getAccumulatedBoxDipole(totalDipole.getArrayPointer());

      curSnapshot->statData[Stats::BOX_DIPOLE_X] = totalDipole(0);
      curSnapshot->statData[Stats::BOX_DIPOLE_Y] = totalDipole(1);
      curSnapshot->statData[Stats::BOX_DIPOLE_Z] = totalDipole(2);
    }

    //store the tau and long range potential
    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
  }


  void ForceManager::postCalculation() {
    SimInfo::MoleculeIterator mi;
    Molecule* mol;
    Molecule::RigidBodyIterator rbIter;
    RigidBody* rb;
    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();

    // collect the atomic forces onto rigid bodies

    for (mol = info_->beginMolecule(mi); mol != NULL;
         mol = info_->nextMolecule(mi)) {
      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
           rb = mol->nextRigidBody(rbIter)) {
        Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
        tau += rbTau;
      }
    }

#ifdef IS_MPI
    Mat3x3d tmpTau(tau);
    MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
                  9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
#endif
    curSnapshot->statData.setTau(tau);
    curSnapshot->statData.setJv(Jv_);
  }

} //end namespace OpenMD
Revision:	1682
Committed:	Tue Feb 28 23:11:22 2012 UTC (13 years, 4 months ago) by chuckv
File size:	12683 byte(s)
Log Message:	Debugging heat flux calculation for rigid bodies.
#	User	Rev	Content
1	gezelter	507	/*
2	gezelter	246	* 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	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	gezelter	246	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	gezelter	246	* 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	gezelter	1390	*
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	chuckv	1671	*
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	gezelter	246	*/
41	chuckv	1671
42	gezelter	507	/**
43			* @file ForceManager.cpp
44			* @author tlin
45			* @date 11/09/2004
46			* @time 10:39am
47			* @version 1.0
48			*/
49	gezelter	246
50			#include "brains/ForceManager.hpp"
51			#include "primitives/Molecule.hpp"
52			#include "UseTheForce/doForces_interface.h"
53	gezelter	1390	#define __OPENMD_C
54	chuckv	664	#include "UseTheForce/DarkSide/fInteractionMap.h"
55	gezelter	246	#include "utils/simError.h"
56	xsun	1215	#include "primitives/Bond.hpp"
57	tim	749	#include "primitives/Bend.hpp"
58	cli2	1275	#include "primitives/Torsion.hpp"
59			#include "primitives/Inversion.hpp"
60	gezelter	1390	namespace OpenMD {
61	gezelter	246
62	gezelter	1464	void ForceManager::calcForces() {
63	chuckv	1671
64	gezelter	246	if (!info_->isFortranInitialized()) {
65	gezelter	507	info_->update();
66	gezelter	246	}
67	chuckv	1671
68	gezelter	246	preCalculation();
69	chuckv	1671
70	gezelter	246	calcShortRangeInteraction();
71
72	gezelter	1464	calcLongRangeInteraction();
73	gezelter	246
74	gezelter	1464	postCalculation();
75	chuckv	1671
76	gezelter	507	}
77	chuckv	1671
78	gezelter	507	void ForceManager::preCalculation() {
79	gezelter	246	SimInfo::MoleculeIterator mi;
80			Molecule* mol;
81			Molecule::AtomIterator ai;
82			Atom* atom;
83			Molecule::RigidBodyIterator rbIter;
84			RigidBody* rb;
85	chuckv	1671
86	gezelter	246	// forces are zeroed here, before any are accumulated.
87			// NOTE: do not rezero the forces in Fortran.
88	chuckv	1671
89			for (mol = info_->beginMolecule(mi); mol != NULL;
90	gezelter	1126	mol = info_->nextMolecule(mi)) {
91	gezelter	507	for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
92	chuckv	1671	atom->zeroForcesAndTorques();
93	gezelter	507	}
94	chuckv	1671
95	gezelter	507	//change the positions of atoms which belong to the rigidbodies
96	chuckv	1671	for (rb = mol->beginRigidBody(rbIter); rb != NULL;
97	gezelter	1126	rb = mol->nextRigidBody(rbIter)) {
98	chuckv	1671	rb->zeroForcesAndTorques();
99			}
100
101	gezelter	246	}
102	chuckv	1671
103	gezelter	1126	// Zero out the stress tensor
104			tau *= 0.0;
105	chuckv	1671
106	gezelter	507	}
107	chuckv	1671
108	gezelter	507	void ForceManager::calcShortRangeInteraction() {
109	gezelter	246	Molecule* mol;
110			RigidBody* rb;
111			Bond* bond;
112			Bend* bend;
113			Torsion* torsion;
114	cli2	1275	Inversion* inversion;
115	gezelter	246	SimInfo::MoleculeIterator mi;
116			Molecule::RigidBodyIterator rbIter;
117			Molecule::BondIterator bondIter;;
118			Molecule::BendIterator bendIter;
119			Molecule::TorsionIterator torsionIter;
120	cli2	1275	Molecule::InversionIterator inversionIter;
121	tim	963	RealType bondPotential = 0.0;
122			RealType bendPotential = 0.0;
123			RealType torsionPotential = 0.0;
124	cli2	1275	RealType inversionPotential = 0.0;
125	gezelter	246
126	chuckv	1671	//calculate short range interactions
127			for (mol = info_->beginMolecule(mi); mol != NULL;
128	gezelter	1126	mol = info_->nextMolecule(mi)) {
129	gezelter	246
130	gezelter	507	//change the positions of atoms which belong to the rigidbodies
131	chuckv	1671	for (rb = mol->beginRigidBody(rbIter); rb != NULL;
132	gezelter	1126	rb = mol->nextRigidBody(rbIter)) {
133			rb->updateAtoms();
134	gezelter	507	}
135	gezelter	246
136	chuckv	1671	for (bond = mol->beginBond(bondIter); bond != NULL;
137	gezelter	1126	bond = mol->nextBond(bondIter)) {
138	tim	749	bond->calcForce();
139			bondPotential += bond->getPotential();
140	gezelter	507	}
141	gezelter	246
142	chuckv	1671	for (bend = mol->beginBend(bendIter); bend != NULL;
143	gezelter	1126	bend = mol->nextBend(bendIter)) {
144	chuckv	1671
145	gezelter	1126	RealType angle;
146			bend->calcForce(angle);
147	chuckv	1671	RealType currBendPot = bend->getPotential();
148
149	gezelter	1126	bendPotential += bend->getPotential();
150			std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
151			if (i == bendDataSets.end()) {
152			BendDataSet dataSet;
153			dataSet.prev.angle = dataSet.curr.angle = angle;
154			dataSet.prev.potential = dataSet.curr.potential = currBendPot;
155			dataSet.deltaV = 0.0;
156			bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
157			}else {
158			i->second.prev.angle = i->second.curr.angle;
159			i->second.prev.potential = i->second.curr.potential;
160			i->second.curr.angle = angle;
161			i->second.curr.potential = currBendPot;
162	chuckv	1671	i->second.deltaV = fabs(i->second.curr.potential -
163	gezelter	1126	i->second.prev.potential);
164			}
165	gezelter	507	}
166	chuckv	1671
167			for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
168	gezelter	1126	torsion = mol->nextTorsion(torsionIter)) {
169	tim	963	RealType angle;
170	gezelter	1126	torsion->calcForce(angle);
171	tim	963	RealType currTorsionPot = torsion->getPotential();
172	gezelter	1126	torsionPotential += torsion->getPotential();
173			std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
174			if (i == torsionDataSets.end()) {
175			TorsionDataSet dataSet;
176			dataSet.prev.angle = dataSet.curr.angle = angle;
177			dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
178			dataSet.deltaV = 0.0;
179			torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
180			}else {
181			i->second.prev.angle = i->second.curr.angle;
182			i->second.prev.potential = i->second.curr.potential;
183			i->second.curr.angle = angle;
184			i->second.curr.potential = currTorsionPot;
185	chuckv	1671	i->second.deltaV = fabs(i->second.curr.potential -
186	gezelter	1126	i->second.prev.potential);
187	chuckv	1671	}
188			}
189	cli2	1275
190	chuckv	1671	for (inversion = mol->beginInversion(inversionIter);
191			inversion != NULL;
192	cli2	1275	inversion = mol->nextInversion(inversionIter)) {
193			RealType angle;
194			inversion->calcForce(angle);
195			RealType currInversionPot = inversion->getPotential();
196			inversionPotential += inversion->getPotential();
197			std::map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
198			if (i == inversionDataSets.end()) {
199			InversionDataSet dataSet;
200			dataSet.prev.angle = dataSet.curr.angle = angle;
201			dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
202			dataSet.deltaV = 0.0;
203			inversionDataSets.insert(std::map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
204			}else {
205			i->second.prev.angle = i->second.curr.angle;
206			i->second.prev.potential = i->second.curr.potential;
207			i->second.curr.angle = angle;
208			i->second.curr.potential = currInversionPot;
209	chuckv	1671	i->second.deltaV = fabs(i->second.curr.potential -
210	cli2	1275	i->second.prev.potential);
211	chuckv	1671	}
212			}
213	gezelter	246	}
214	chuckv	1671
215			RealType shortRangePotential = bondPotential + bendPotential +
216			torsionPotential + inversionPotential;
217	gezelter	246	Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
218			curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
219	tim	665	curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
220			curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
221			curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
222	cli2	1275	curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
223	chuckv	1671
224	gezelter	507	}
225	chuckv	1671
226	gezelter	1464	void ForceManager::calcLongRangeInteraction() {
227	gezelter	246	Snapshot* curSnapshot;
228			DataStorage* config;
229	tim	963	RealType* frc;
230			RealType* pos;
231	chuckv	1671	RealType* vel;
232	tim	963	RealType* trq;
233			RealType* A;
234			RealType* electroFrame;
235			RealType* rc;
236	chuckv	1682	RealType* vc;
237	chuckv	1245	RealType* particlePot;
238	chuckv	1671
239	gezelter	246	//get current snapshot from SimInfo
240			curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
241	chuckv	1671
242	gezelter	246	//get array pointers
243			config = &(curSnapshot->atomData);
244			frc = config->getArrayPointer(DataStorage::dslForce);
245			pos = config->getArrayPointer(DataStorage::dslPosition);
246	chuckv	1671	vel = config->getArrayPointer(DataStorage::dslVelocity);
247	gezelter	246	trq = config->getArrayPointer(DataStorage::dslTorque);
248			A = config->getArrayPointer(DataStorage::dslAmat);
249			electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
250	chuckv	1245	particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
251	gezelter	246
252			//calculate the center of mass of cutoff group
253			SimInfo::MoleculeIterator mi;
254			Molecule* mol;
255			Molecule::CutoffGroupIterator ci;
256			CutoffGroup* cg;
257			Vector3d com;
258	chuckv	1682	Vector3d comv;
259	gezelter	246	std::vector<Vector3d> rcGroup;
260	chuckv	1682	std::vector<Vector3d> vcGroup;
261	chuckv	1671
262	gezelter	246	if(info_->getNCutoffGroups() > 0){
263	chuckv	1671
264			for (mol = info_->beginMolecule(mi); mol != NULL;
265	gezelter	1126	mol = info_->nextMolecule(mi)) {
266	chuckv	1671	for(cg = mol->beginCutoffGroup(ci); cg != NULL;
267	gezelter	1126	cg = mol->nextCutoffGroup(ci)) {
268	chuckv	1682	cg->getCOM(com);
269			cg->getCOMV(comv);
270			rcGroup.push_back(com);
271			vcGroup.push_back(comv);
272	gezelter	246	}
273	gezelter	507	}// end for (mol)
274	chuckv	1671
275	gezelter	507	rc = rcGroup[0].getArrayPointer();
276	chuckv	1682	vc = vcGroup[0].getArrayPointer();
277	gezelter	246	} else {
278	chuckv	1671	// center of mass of the group is the same as position of the atom
279	gezelter	1126	// if cutoff group does not exist
280	gezelter	507	rc = pos;
281	chuckv	1682	vc = vel;
282	gezelter	246	}
283	chuckv	1671
284	gezelter	246	//initialize data before passing to fortran
285	tim	963	RealType longRangePotential[LR_POT_TYPES];
286			RealType lrPot = 0.0;
287	chrisfen	998	Vector3d totalDipole;
288	chuckv	1681	Jv_ = 0.0;
289	chuckv	1671
290
291
292	gezelter	246	int isError = 0;
293
294	chuckv	664	for (int i=0; i<LR_POT_TYPES;i++){
295			longRangePotential[i]=0.0; //Initialize array
296			}
297	chuckv	1671
298	xsun	1215	doForceLoop(pos,
299	chuckv	1671	vel,
300	xsun	1215	rc,
301	chuckv	1682	vc,
302	xsun	1215	A,
303			electroFrame,
304			frc,
305			trq,
306	chuckv	1671	tau.getArrayPointer(),
307			Jv_.getArrayPointer(),
308			longRangePotential,
309	chuckv	1245	particlePot,
310	xsun	1215	&isError );
311	chuckv	1671
312	gezelter	246	if( isError ){
313	gezelter	507	sprintf( painCave.errMsg,
314	chuckv	1671	"Error returned from the fortran force calculation.\n" );
315	gezelter	507	painCave.isFatal = 1;
316			simError();
317	gezelter	246	}
318	chuckv	664	for (int i=0; i<LR_POT_TYPES;i++){
319			lrPot += longRangePotential[i]; //Quick hack
320			}
321	chuckv	1671
322	chrisfen	998	// grab the simulation box dipole moment if specified
323			if (info_->getCalcBoxDipole()){
324			getAccumulatedBoxDipole(totalDipole.getArrayPointer());
325	chuckv	1671
326	chrisfen	998	curSnapshot->statData[Stats::BOX_DIPOLE_X] = totalDipole(0);
327			curSnapshot->statData[Stats::BOX_DIPOLE_Y] = totalDipole(1);
328			curSnapshot->statData[Stats::BOX_DIPOLE_Z] = totalDipole(2);
329			}
330	chuckv	1671
331			//store the tau and long range potential
332	chuckv	664	curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
333	chrisfen	691	curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
334	tim	681	curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
335	gezelter	507	}
336	gezelter	246
337	chuckv	1671
338	gezelter	1464	void ForceManager::postCalculation() {
339	gezelter	246	SimInfo::MoleculeIterator mi;
340			Molecule* mol;
341			Molecule::RigidBodyIterator rbIter;
342			RigidBody* rb;
343	gezelter	1126	Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
344	chuckv	1671
345	gezelter	246	// collect the atomic forces onto rigid bodies
346	chuckv	1671
347			for (mol = info_->beginMolecule(mi); mol != NULL;
348	gezelter	1126	mol = info_->nextMolecule(mi)) {
349	chuckv	1671	for (rb = mol->beginRigidBody(rbIter); rb != NULL;
350			rb = mol->nextRigidBody(rbIter)) {
351	gezelter	1464	Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
352			tau += rbTau;
353	gezelter	507	}
354	gezelter	1126	}
355	chuckv	1671
356	gezelter	1126	#ifdef IS_MPI
357	gezelter	1464	Mat3x3d tmpTau(tau);
358	chuckv	1671	MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
359	gezelter	1464	9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
360	gezelter	1126	#endif
361	gezelter	1464	curSnapshot->statData.setTau(tau);
362	chuckv	1671	curSnapshot->statData.setJv(Jv_);
363	gezelter	507	}
364	gezelter	246
365	gezelter	1390	} //end namespace OpenMD