9 |
|
#include "SRI.hpp" |
10 |
|
#include "simError.h" |
11 |
|
|
12 |
+ |
#include <fortranWrappers.hpp> |
13 |
+ |
|
14 |
|
#ifdef IS_MPI |
15 |
|
#include "mpiForceField.h" |
16 |
|
#endif // is_mpi |
86 |
|
|
87 |
|
} // namespace |
88 |
|
|
87 |
– |
|
88 |
– |
// declaration of functions needed to wrap the fortran module |
89 |
– |
|
90 |
– |
|
91 |
– |
|
92 |
– |
TraPPE_ExFF* currentTPEwrap; |
93 |
– |
|
89 |
|
using namespace TPE; |
90 |
|
|
91 |
|
|
103 |
|
char errMsg[1000]; |
104 |
|
|
105 |
|
// do the funtion wrapping |
106 |
< |
currentTPEwrap = this; |
112 |
< |
wrapMe(); |
106 |
> |
wrapMeFF( this ); |
107 |
|
|
108 |
|
|
109 |
|
#ifdef IS_MPI |
217 |
|
|
218 |
|
ffPath = getenv( ffPath_env ); |
219 |
|
if( ffPath == NULL ) { |
220 |
< |
sprintf( painCave.errMsg, |
227 |
< |
"Error opening the force field parameter file: %s\n" |
228 |
< |
"Have you tried setting the FORCE_PARAM_PATH environment " |
229 |
< |
"vairable?\n", |
230 |
< |
fileName ); |
231 |
< |
painCave.isFatal = 1; |
232 |
< |
simError(); |
220 |
> |
STR_DEFINE(ffPath, FRC_PATH ); |
221 |
|
} |
222 |
|
|
223 |
|
|
263 |
|
#endif // is_mpi |
264 |
|
} |
265 |
|
|
266 |
+ |
void TraPPE_ExFF::doForces( int calcPot ){ |
267 |
|
|
268 |
< |
void TraPPE_ExFF::wrapMe( void ){ |
269 |
< |
|
270 |
< |
char* currentFF = "TraPPE_Ex"; |
271 |
< |
int isError = 0; |
272 |
< |
|
284 |
< |
forcefactory_( currentFF, &isError, TPEfunctionWrapper, strlen(currentFF) ); |
268 |
> |
int i, isError; |
269 |
> |
double* frc; |
270 |
> |
double* pos; |
271 |
> |
double* tau; |
272 |
> |
short int passedCalcPot = (short int)calcPot; |
273 |
|
|
274 |
< |
if( isError ){ |
275 |
< |
|
274 |
> |
// forces are zeroed here, before any are acumulated. |
275 |
> |
// NOTE: do not rezero the forces in Fortran. |
276 |
> |
|
277 |
> |
for(i=0; i<entry_plug->n_atoms; i++){ |
278 |
> |
entry_plug->atoms[i]->zeroForces(); |
279 |
> |
} |
280 |
> |
|
281 |
> |
frc = Atom::getFrcArray(); |
282 |
> |
pos = Atom::getPosArray(); |
283 |
> |
tau = entry_plug->tau; |
284 |
> |
|
285 |
> |
isError = 0; |
286 |
> |
fortranForceLoop( pos, frc, &(entry_plug->lrPot), tau, |
287 |
> |
&passedCalcPot, &isError ); |
288 |
> |
|
289 |
> |
|
290 |
> |
if( isError ){ |
291 |
|
sprintf( painCave.errMsg, |
292 |
< |
"TraPPE_ExFF error: an error was returned from fortran when the " |
290 |
< |
"the functions were being wrapped.\n" ); |
292 |
> |
"Error returned from the fortran force calculation.\n" ); |
293 |
|
painCave.isFatal = 1; |
294 |
|
simError(); |
295 |
|
} |
296 |
|
|
297 |
|
#ifdef IS_MPI |
298 |
< |
sprintf( checkPointMsg, "TraPPE_ExFF functions succesfully wrapped." ); |
298 |
> |
sprintf( checkPointMsg, |
299 |
> |
"successfully returned from the force calculation.\n" ); |
300 |
|
MPIcheckPoint(); |
301 |
|
#endif // is_mpi |
302 |
+ |
|
303 |
|
} |
304 |
|
|
305 |
+ |
void TraPPE_ExFF::initFortran( void ){ |
306 |
+ |
|
307 |
+ |
int nLocal = entry_plug->n_atoms; |
308 |
+ |
int *ident; |
309 |
+ |
int isError; |
310 |
+ |
int i; |
311 |
|
|
312 |
< |
void TPEfunctionWrapper( void (*p1)( int* ident, double* mass, double* epslon, |
313 |
< |
double* sigma, int* isDipole, |
314 |
< |
int* isSSD, double* dipole, double* w0, |
315 |
< |
double* v0, int* isError ), |
316 |
< |
void (*p2)( int*, int*, int* ), |
317 |
< |
void (*p3)( double*,double*,double*,double*, |
318 |
< |
short int*, int* ) ){ |
312 |
> |
ident = new int[nLocal]; |
313 |
> |
|
314 |
> |
for(i=0; i<nLocal; i++){ |
315 |
> |
ident[i] = entry_plug->atoms[i]->getIdent(); |
316 |
> |
} |
317 |
> |
|
318 |
> |
isError = 0; |
319 |
> |
initfortran( &nLocal, ident, &isError ); |
320 |
|
|
321 |
< |
|
322 |
< |
newTPEtype = p1; |
323 |
< |
initTPEfortran = p2; |
324 |
< |
currentTPEwrap->setTPEfortran( p3 ); |
321 |
> |
if(isError){ |
322 |
> |
sprintf( painCave.errMsg, |
323 |
> |
"TraPPE_ExFF error: There was an error initializing the component list in fortran.\n" ); |
324 |
> |
painCave.isFatal = 1; |
325 |
> |
simError(); |
326 |
> |
} |
327 |
> |
|
328 |
> |
|
329 |
> |
#ifdef IS_MPI |
330 |
> |
sprintf( checkPointMsg, "TraPPE_ExFF successfully initialized the fortran component list.\n" ); |
331 |
> |
MPIcheckPoint(); |
332 |
> |
#endif // is_mpi |
333 |
> |
|
334 |
> |
delete[] ident; |
335 |
> |
|
336 |
|
} |
337 |
|
|
338 |
|
|
339 |
|
|
340 |
+ |
|
341 |
|
void TraPPE_ExFF::initializeAtoms( void ){ |
342 |
|
|
343 |
|
class LinkedType { |
563 |
|
// call new A_types in fortran |
564 |
|
|
565 |
|
int isError; |
566 |
+ |
|
567 |
+ |
// dummy variables |
568 |
+ |
|
569 |
+ |
int isGB = 0; |
570 |
+ |
int isLJ = 1; |
571 |
+ |
|
572 |
+ |
|
573 |
|
currentAtomType = headAtomType; |
574 |
|
while( currentAtomType != NULL ){ |
575 |
|
|
579 |
|
&(currentAtomType->mass), |
580 |
|
&(currentAtomType->epslon), |
581 |
|
&(currentAtomType->sigma), |
582 |
< |
&(currentAtomType->isDipole), |
582 |
> |
&isLJ, |
583 |
|
&(currentAtomType->isSSD), |
584 |
< |
&(currentAtomType->dipole), |
584 |
> |
&(currentAtomType->isDipole), |
585 |
> |
&isGB, |
586 |
|
&(currentAtomType->w0), |
587 |
|
&(currentAtomType->v0), |
588 |
+ |
&(currentAtomType->dipole), |
589 |
|
&isError ); |
590 |
|
if( isError ){ |
591 |
|
sprintf( painCave.errMsg, |
707 |
|
MPIcheckPoint(); |
708 |
|
#endif // is_mpi |
709 |
|
|
710 |
+ |
initFortran(); |
711 |
+ |
entry_plug->refreshSim(); |
712 |
+ |
|
713 |
|
} |
714 |
|
|
715 |
|
void TraPPE_ExFF::initializeBonds( bond_pair* the_bonds ){ |
1003 |
|
// if things go well, last will be set to 0 |
1004 |
|
|
1005 |
|
QuadraticBend* qBend; |
1006 |
+ |
GhostBend* gBend; |
1007 |
|
SRI **the_sris; |
1008 |
|
Atom** the_atoms; |
1009 |
|
int nBends; |
1104 |
|
|
1105 |
|
atomA = the_atoms[a]->getType(); |
1106 |
|
atomB = the_atoms[b]->getType(); |
1107 |
< |
atomC = the_atoms[c]->getType(); |
1107 |
> |
|
1108 |
> |
if( the_bends[i].isGhost ) atomC = "GHOST"; |
1109 |
> |
else atomC = the_atoms[c]->getType(); |
1110 |
> |
|
1111 |
|
currentBendType = headBendType->find( atomA, atomB, atomC ); |
1112 |
|
if( currentBendType == NULL ){ |
1113 |
|
sprintf( painCave.errMsg, "BendType error, %s - %s - %s not found" |
1120 |
|
if( !strcmp( currentBendType->type, "quadratic" ) ){ |
1121 |
|
|
1122 |
|
index = i + entry_plug->n_bonds; |
1123 |
< |
qBend = new QuadraticBend( *the_atoms[a], |
1124 |
< |
*the_atoms[b], |
1125 |
< |
*the_atoms[c] ); |
1126 |
< |
qBend->setConstants( currentBendType->k1, |
1127 |
< |
currentBendType->k2, |
1128 |
< |
currentBendType->k3, |
1129 |
< |
currentBendType->t0 ); |
1130 |
< |
the_sris[index] = qBend; |
1123 |
> |
|
1124 |
> |
if( the_bends[i].isGhost){ |
1125 |
> |
|
1126 |
> |
if( the_bends[i].ghost == b ){ |
1127 |
> |
// do nothing |
1128 |
> |
} |
1129 |
> |
else if( the_bends[i].ghost == a ){ |
1130 |
> |
c = a; |
1131 |
> |
a = b; |
1132 |
> |
b = a; |
1133 |
> |
} |
1134 |
> |
else{ |
1135 |
> |
sprintf( painCave.errMsg, |
1136 |
> |
"BendType error, %s - %s - %s,\n" |
1137 |
> |
" --> central atom is not " |
1138 |
> |
"correctly identified with the " |
1139 |
> |
"\"ghostVectorSource = \" tag.\n", |
1140 |
> |
atomA, atomB, atomC ); |
1141 |
> |
painCave.isFatal = 1; |
1142 |
> |
simError(); |
1143 |
> |
} |
1144 |
> |
|
1145 |
> |
gBend = new GhostBend( *the_atoms[a], |
1146 |
> |
*the_atoms[b] ); |
1147 |
> |
gBend->setConstants( currentBendType->k1, |
1148 |
> |
currentBendType->k2, |
1149 |
> |
currentBendType->k3, |
1150 |
> |
currentBendType->t0 ); |
1151 |
> |
the_sris[index] = gBend; |
1152 |
> |
} |
1153 |
> |
else{ |
1154 |
> |
qBend = new QuadraticBend( *the_atoms[a], |
1155 |
> |
*the_atoms[b], |
1156 |
> |
*the_atoms[c] ); |
1157 |
> |
qBend->setConstants( currentBendType->k1, |
1158 |
> |
currentBendType->k2, |
1159 |
> |
currentBendType->k3, |
1160 |
> |
currentBendType->t0 ); |
1161 |
> |
the_sris[index] = qBend; |
1162 |
> |
} |
1163 |
|
} |
1164 |
|
} |
1165 |
|
|