src/brains/ForceManager.cpp

/*
 * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
 * 
 * Contact: oopse@oopse.org
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */
 
 /**
  * @file ForceManager.cpp
  * @author tlin
  * @date 11/09/2004
  * @time 10:39am
  * @version 1.0
  */

#include "brains/ForceManager.hpp"

namespace oopse {

void ForceManager::calcForces(bool needPotential, bool needStress) {

    Snapshot* curSnapshot;
    DataStorage* config;
    double* frc;
    double* pos;
    double* trq;
    double* A;
    double* u_l;
    double* rc;

    typename SimInfo::MoleculeIterator mi;
    Molecule* mol;
    std::vector<Atom*>::iterator ai;
    Atom* atom;
    std::vector<Atom*>::iterator ci;
    CutoffGroup* cg;
    Vector3d com;
    std::vector<Vector3d> rcGroup;

    double longRangePotential;
    double tau[9];
    short int passedCalcPot;
    short int passedCalcStress;
    int isError;     

    // 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->zeroForces();
        }
    }// end for (mol)

#ifdef PROFILE
    startProfile(pro7);
#endif
    //calculate short range interactions    
    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
        mol->calcForces();
    }// end for (mol)

#ifdef PROFILE
    endProfile( pro7 );
#endif

    //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);
    trq = config->getArrayPointer(DataStorage::dslTorque);
    A   = config->getArrayPointer(DataStorage::dslAmat);
    u_l = config->getArrayPointer(DataStorage::dslUnitVector);

    //calculate the center of mass of cutoff group
    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)
            rcGroup.push_back(com);
        }
    }// end for (mol)
       
        rc = rcGroup[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;
    }
  

    //initialize data before passing to fortran
    longRangePotential = 0.0;
    for (int i=0; i<9; i++) {
        tau[i] = 0.0;
    }    
    passedCalcPot = (short int)needPotential;
    passedCalcStress = (short int)needStress;
    isError = 0;

#ifdef PROFILE
    startProfile(pro8);
#endif

    doForceLoop( pos,
            rc,
            A,
            u_l,
            frc,
            trq,
            tau,
            &longRangePotential, 
            &passedCalcPot,
            &passedCalcStress,
            &isError );

#ifdef PROFILE
    endProfile(pro8);
#endif

    if( isError ){
        sprintf( painCave.errMsg,
             "Error returned from the fortran force calculation.\n" );
        painCave.isFatal = 1;
        simError();
    }

    //store the tau and long range potential    
    curSnapshot->statData[Stats::LONGRANGEPOTENTIAL] = longRangePotential;
    curSnapshot->statData[Stats::TAUXX] = tau[0];
    curSnapshot->statData[Stats::TAUXY] = tau[1];
    curSnapshot->statData[Stats::TAUXZ] = tau[2];
    curSnapshot->statData[Stats::TAUYX] = tau[3];
    curSnapshot->statData[Stats::TAUYY] = tau[4];
    curSnapshot->statData[Stats::TAUYZ] = tau[5];
    curSnapshot->statData[Stats::TAUZX] = tau[6];
    curSnapshot->statData[Stats::TAUZY] = tau[7];
    curSnapshot->statData[Stats::TAUZZ] = tau[8];

    // collect the atomic forces onto rigid bodies
    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
        mol->atoms2rigidBodies();
    }

#ifdef IS_MPI
    sprintf( checkPointMsg,
       "returned from the force calculation.\n" );
    MPIcheckPoint();
#endif // is_mpi

}


} //end namespace oopse
Revision:	1726
Committed:	Wed Nov 10 22:50:03 2004 UTC (19 years, 8 months ago) by tim
File size:	5228 byte(s)
Log Message:	more work on SimInfo
#	Content
1	/*
2	* Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project
3	*
4	* Contact: oopse@oopse.org
5	*
6	* This program is free software; you can redistribute it and/or
7	* modify it under the terms of the GNU Lesser General Public License
8	* as published by the Free Software Foundation; either version 2.1
9	* of the License, or (at your option) any later version.
10	* All we ask is that proper credit is given for our work, which includes
11	* - but is not limited to - adding the above copyright notice to the beginning
12	* of your source code files, and to any copyright notice that you may distribute
13	* with programs based on this work.
14	*
15	* This program is distributed in the hope that it will be useful,
16	* but WITHOUT ANY WARRANTY; without even the implied warranty of
17	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18	* GNU Lesser General Public License for more details.
19	*
20	* You should have received a copy of the GNU Lesser General Public License
21	* along with this program; if not, write to the Free Software
22	* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
23	*
24	*/
25
26	/**
27	* @file ForceManager.cpp
28	* @author tlin
29	* @date 11/09/2004
30	* @time 10:39am
31	* @version 1.0
32	*/
33
34	#include "brains/ForceManager.hpp"
35
36	namespace oopse {
37
38	void ForceManager::calcForces(bool needPotential, bool needStress) {
39
40	Snapshot* curSnapshot;
41	DataStorage* config;
42	double* frc;
43	double* pos;
44	double* trq;
45	double* A;
46	double* u_l;
47	double* rc;
48
49	typename SimInfo::MoleculeIterator mi;
50	Molecule* mol;
51	std::vector<Atom*>::iterator ai;
52	Atom* atom;
53	std::vector<Atom*>::iterator ci;
54	CutoffGroup* cg;
55	Vector3d com;
56	std::vector<Vector3d> rcGroup;
57
58	double longRangePotential;
59	double tau[9];
60	short int passedCalcPot;
61	short int passedCalcStress;
62	int isError;
63
64	// forces are zeroed here, before any are accumulated.
65	// NOTE: do not rezero the forces in Fortran.
66	for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
67	for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
68	atom->zeroForces();
69	}
70	}// end for (mol)
71
72	#ifdef PROFILE
73	startProfile(pro7);
74	#endif
75	//calculate short range interactions
76	for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
77	mol->calcForces();
78	}// end for (mol)
79
80	#ifdef PROFILE
81	endProfile( pro7 );
82	#endif
83
84	//get current snapshot from SimInfo
85	curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot()
86
87	//get array pointers
88	config = &(curSnapshot->atomData);
89	frc = config->getArrayPointer(DataStorage::dslForce);
90	pos = config->getArrayPointer(DataStorage::dslPosition);
91	trq = config->getArrayPointer(DataStorage::dslTorque);
92	A = config->getArrayPointer(DataStorage::dslAmat);
93	u_l = config->getArrayPointer(DataStorage::dslUnitVector);
94
95	//calculate the center of mass of cutoff group
96	if(info_->getNCutoffGroups() > 0){
97
98	for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
99	for(cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
100	cg->getCOM(com)
101	rcGroup.push_back(com);
102	}
103	}// end for (mol)
104
105	rc = rcGroup[0]->getArrayPointer();
106	}
107	else{
108	// center of mass of the group is the same as position of the atom if cutoff group does not exist
109	rc = pos;
110	}
111
112
113	//initialize data before passing to fortran
114	longRangePotential = 0.0;
115	for (int i=0; i<9; i++) {
116	tau[i] = 0.0;
117	}
118	passedCalcPot = (short int)needPotential;
119	passedCalcStress = (short int)needStress;
120	isError = 0;
121
122	#ifdef PROFILE
123	startProfile(pro8);
124	#endif
125
126	doForceLoop( pos,
127	rc,
128	A,
129	u_l,
130	frc,
131	trq,
132	tau,
133	&longRangePotential,
134	&passedCalcPot,
135	&passedCalcStress,
136	&isError );
137
138	#ifdef PROFILE
139	endProfile(pro8);
140	#endif
141
142	if( isError ){
143	sprintf( painCave.errMsg,
144	"Error returned from the fortran force calculation.\n" );
145	painCave.isFatal = 1;
146	simError();
147	}
148
149	//store the tau and long range potential
150	curSnapshot->statData[Stats::LONGRANGEPOTENTIAL] = longRangePotential;
151	curSnapshot->statData[Stats::TAUXX] = tau[0];
152	curSnapshot->statData[Stats::TAUXY] = tau[1];
153	curSnapshot->statData[Stats::TAUXZ] = tau[2];
154	curSnapshot->statData[Stats::TAUYX] = tau[3];
155	curSnapshot->statData[Stats::TAUYY] = tau[4];
156	curSnapshot->statData[Stats::TAUYZ] = tau[5];
157	curSnapshot->statData[Stats::TAUZX] = tau[6];
158	curSnapshot->statData[Stats::TAUZY] = tau[7];
159	curSnapshot->statData[Stats::TAUZZ] = tau[8];
160
161	// collect the atomic forces onto rigid bodies
162	for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
163	mol->atoms2rigidBodies();
164	}
165
166	#ifdef IS_MPI
167	sprintf( checkPointMsg,
168	"returned from the force calculation.\n" );
169	MPIcheckPoint();
170	#endif // is_mpi
171
172	}
173
174
175	} //end namespace oopse