73 |
|
|
74 |
|
void ForceFields::doForces( int calcPot, int calcStress ){ |
75 |
|
|
76 |
< |
int i, isError; |
76 |
> |
int i, j, isError; |
77 |
|
double* frc; |
78 |
|
double* pos; |
79 |
|
double* trq; |
81 |
|
double* u_l; |
82 |
|
double* rc; |
83 |
|
double* massRatio; |
84 |
+ |
double factor; |
85 |
|
SimState* config; |
86 |
|
|
87 |
+ |
Molecule* myMols; |
88 |
+ |
Atom** myAtoms; |
89 |
+ |
int numAtom; |
90 |
+ |
int curIndex; |
91 |
+ |
double mtot; |
92 |
+ |
int numMol; |
93 |
+ |
int numCutoffGroups; |
94 |
+ |
CutoffGroup* myCutoffGroup; |
95 |
+ |
vector<CutoffGroup*>::iterator iterCutoff; |
96 |
+ |
double com[3]; |
97 |
+ |
vector<double> rcGroup; |
98 |
+ |
|
99 |
|
short int passedCalcPot = (short int)calcPot; |
100 |
|
short int passedCalcStress = (short int)calcStress; |
101 |
|
|
113 |
|
for(i=0; i<entry_plug->n_mol; i++ ){ |
114 |
|
// CalcForces in molecules takes care of mapping rigid body coordinates |
115 |
|
// into atomic coordinates |
116 |
< |
entry_plug->molecules[i].calcForces(); |
116 |
> |
entry_plug->molecules[i].calcForces(); |
117 |
|
} |
118 |
|
|
119 |
|
#ifdef PROFILE |
127 |
|
trq = config->getTrqArray(); |
128 |
|
A = config->getAmatArray(); |
129 |
|
u_l = config->getUlArray(); |
117 |
– |
rc = config->getRcArray(); |
118 |
– |
massRatio = config->getMassRatioArray(); |
130 |
|
|
131 |
+ |
if(entry_plug->haveCutoffGroups){ |
132 |
+ |
myMols = entry_plug->molecules; |
133 |
+ |
numMol = entry_plug->n_mol; |
134 |
+ |
for(int i = 0; i < numMol; i++){ |
135 |
+ |
|
136 |
+ |
numCutoffGroups = myMols[i].getNCutoffGroups(); |
137 |
+ |
for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); myCutoffGroup != NULL; |
138 |
+ |
myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){ |
139 |
+ |
//get center of mass of the cutoff group |
140 |
+ |
myCutoffGroup->getCOM(com); |
141 |
+ |
|
142 |
+ |
rcGroup.push_back(com[0]); |
143 |
+ |
rcGroup.push_back(com[1]); |
144 |
+ |
rcGroup.push_back(com[2]); |
145 |
+ |
|
146 |
+ |
}// end for(myCutoffGroup) |
147 |
+ |
|
148 |
+ |
}//end for(int i = 0) |
149 |
+ |
|
150 |
+ |
rc = &rcGroup[0]; |
151 |
+ |
} |
152 |
+ |
else{ |
153 |
+ |
// center of mass of the group is the same as position of the atom if cutoff group does not exist |
154 |
+ |
rc = pos; |
155 |
+ |
} |
156 |
+ |
|
157 |
+ |
|
158 |
+ |
|
159 |
|
isError = 0; |
160 |
|
entry_plug->lrPot = 0.0; |
161 |
|
|
198 |
|
} |
199 |
|
|
200 |
|
|
201 |
+ |
if (entry_plug->useThermInt) { |
202 |
+ |
|
203 |
+ |
factor = pow(entry_plug->thermIntLambda, entry_plug->thermIntK); |
204 |
+ |
for (i=0; i < entry_plug->integrableObjects.size(); i++) { |
205 |
+ |
for (j=0; j< 3; j++) |
206 |
+ |
frc[3*i + j] *= factor; |
207 |
+ |
if (entry_plug->integrableObjects[i]->isDirectional()) { |
208 |
+ |
for (j=0; j< 3; j++) |
209 |
+ |
trq[3*i + j] *= factor; |
210 |
+ |
} |
211 |
+ |
} |
212 |
+ |
entry_plug->vRaw = entry_plug->lrPot; |
213 |
+ |
entry_plug->lrPot *= factor; |
214 |
+ |
entry_plug->lrPot += entry_plug->restraint->Calc_Restraint_Forces(entry_plug->integrableObjects); |
215 |
+ |
entry_plug->vHarm = entry_plug->restraint->getVharm(); |
216 |
+ |
} |
217 |
+ |
|
218 |
|
#ifdef IS_MPI |
219 |
|
sprintf( checkPointMsg, |
220 |
|
"returned from the force calculation.\n" ); |
246 |
|
#endif // is_mpi |
247 |
|
|
248 |
|
} |
249 |
+ |
|
250 |
+ |
|
251 |
+ |
void ForceFields::initRestraints(){ |
252 |
+ |
int i; |
253 |
+ |
// store the initial info. |
254 |
+ |
// set the omega values to zero |
255 |
+ |
for (i=0; i<entry_plug->integrableObjects.size(); i++) |
256 |
+ |
entry_plug->integrableObjects[i]->setZangle( 0.0 ); |
257 |
+ |
|
258 |
+ |
entry_plug->restraint->Store_Init_Info(entry_plug->integrableObjects); |
259 |
+ |
|
260 |
+ |
} |
261 |
+ |
|
262 |
+ |
void ForceFields::dumpzAngle(){ |
263 |
+ |
|
264 |
+ |
// store the initial info. |
265 |
+ |
entry_plug->restraint->Write_zAngle_File(entry_plug->integrableObjects); |
266 |
+ |
|
267 |
+ |
} |