| 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 |
| 269 |
|
#endif // is_mpi |
| 270 |
|
} |
| 271 |
|
|
| 272 |
+ |
void TraPPE_ExFF::doForces( int calcPot ){ |
| 273 |
|
|
| 274 |
< |
void TraPPE_ExFF::wrapMe( void ){ |
| 275 |
< |
|
| 276 |
< |
char* currentFF = "TraPPE_Ex"; |
| 277 |
< |
int isError = 0; |
| 278 |
< |
|
| 284 |
< |
forcefactory_( currentFF, &isError, TPEfunctionWrapper, strlen(currentFF) ); |
| 274 |
> |
int i, isError; |
| 275 |
> |
double* frc; |
| 276 |
> |
double* pos; |
| 277 |
> |
double* tau; |
| 278 |
> |
short int passedCalcPot = (short int)calcPot; |
| 279 |
|
|
| 280 |
< |
if( isError ){ |
| 281 |
< |
|
| 280 |
> |
// forces are zeroed here, before any are acumulated. |
| 281 |
> |
// NOTE: do not rezero the forces in Fortran. |
| 282 |
> |
|
| 283 |
> |
for(i=0; i<entry_plug->n_atoms; i++){ |
| 284 |
> |
entry_plug->atoms[i]->zeroForces(); |
| 285 |
> |
} |
| 286 |
> |
|
| 287 |
> |
frc = Atom::getFrcArray(); |
| 288 |
> |
pos = Atom::getPosArray(); |
| 289 |
> |
tau = entry_plug->tau; |
| 290 |
> |
|
| 291 |
> |
isError = 0; |
| 292 |
> |
fortranForceLoop( pos, frc, &(entry_plug->lrPot), tau, |
| 293 |
> |
&passedCalcPot, &isError ); |
| 294 |
> |
|
| 295 |
> |
|
| 296 |
> |
if( isError ){ |
| 297 |
|
sprintf( painCave.errMsg, |
| 298 |
< |
"TraPPE_ExFF error: an error was returned from fortran when the " |
| 290 |
< |
"the functions were being wrapped.\n" ); |
| 298 |
> |
"Error returned from the fortran force calculation.\n" ); |
| 299 |
|
painCave.isFatal = 1; |
| 300 |
|
simError(); |
| 301 |
|
} |
| 302 |
|
|
| 303 |
|
#ifdef IS_MPI |
| 304 |
< |
sprintf( checkPointMsg, "TraPPE_ExFF functions succesfully wrapped." ); |
| 304 |
> |
sprintf( checkPointMsg, |
| 305 |
> |
"successfully returned from the force calculation.\n" ); |
| 306 |
|
MPIcheckPoint(); |
| 307 |
|
#endif // is_mpi |
| 308 |
+ |
|
| 309 |
|
} |
| 310 |
|
|
| 311 |
+ |
void TraPPE_ExFF::initFortran( void ){ |
| 312 |
+ |
|
| 313 |
+ |
int nLocal = entry_plug->n_atoms; |
| 314 |
+ |
int *ident; |
| 315 |
+ |
int isError; |
| 316 |
+ |
int i; |
| 317 |
|
|
| 318 |
< |
void TPEfunctionWrapper( void (*p1)( int* ident, double* mass, double* epslon, |
| 319 |
< |
double* sigma, int* isDipole, |
| 320 |
< |
int* isSSD, double* dipole, double* w0, |
| 321 |
< |
double* v0, int* isError ), |
| 322 |
< |
void (*p2)( int*, int*, int* ), |
| 323 |
< |
void (*p3)( double*,double*,double*,double*, |
| 324 |
< |
short int*, int* ) ){ |
| 318 |
> |
ident = new int[nLocal]; |
| 319 |
> |
|
| 320 |
> |
for(i=0; i<nLocal; i++){ |
| 321 |
> |
ident[i] = entry_plug->atoms[i]->getIdent(); |
| 322 |
> |
} |
| 323 |
> |
|
| 324 |
> |
isError = 0; |
| 325 |
> |
initfortran( &nLocal, ident, &isError ); |
| 326 |
|
|
| 327 |
+ |
if(isError){ |
| 328 |
+ |
sprintf( painCave.errMsg, |
| 329 |
+ |
"TraPPE_ExFF error: There was an error initializing the component list in fortran.\n" ); |
| 330 |
+ |
painCave.isFatal = 1; |
| 331 |
+ |
simError(); |
| 332 |
+ |
} |
| 333 |
+ |
|
| 334 |
|
|
| 335 |
< |
newTPEtype = p1; |
| 336 |
< |
initTPEfortran = p2; |
| 337 |
< |
currentTPEwrap->setTPEfortran( p3 ); |
| 335 |
> |
#ifdef IS_MPI |
| 336 |
> |
sprintf( checkPointMsg, "TraPPE_ExFF successfully initialized the fortran component list.\n" ); |
| 337 |
> |
MPIcheckPoint(); |
| 338 |
> |
#endif // is_mpi |
| 339 |
> |
|
| 340 |
> |
delete[] ident; |
| 341 |
> |
|
| 342 |
|
} |
| 343 |
|
|
| 344 |
|
|
| 345 |
|
|
| 346 |
+ |
|
| 347 |
|
void TraPPE_ExFF::initializeAtoms( void ){ |
| 348 |
|
|
| 349 |
|
class LinkedType { |
| 569 |
|
// call new A_types in fortran |
| 570 |
|
|
| 571 |
|
int isError; |
| 572 |
+ |
|
| 573 |
+ |
// dummy variables |
| 574 |
+ |
|
| 575 |
+ |
int isGB = 0; |
| 576 |
+ |
int isLJ = 1; |
| 577 |
+ |
|
| 578 |
+ |
|
| 579 |
|
currentAtomType = headAtomType; |
| 580 |
|
while( currentAtomType != NULL ){ |
| 581 |
|
|
| 585 |
|
&(currentAtomType->mass), |
| 586 |
|
&(currentAtomType->epslon), |
| 587 |
|
&(currentAtomType->sigma), |
| 588 |
< |
&(currentAtomType->isDipole), |
| 588 |
> |
&isLJ, |
| 589 |
|
&(currentAtomType->isSSD), |
| 590 |
< |
&(currentAtomType->dipole), |
| 590 |
> |
&(currentAtomType->isDipole), |
| 591 |
> |
&isGB, |
| 592 |
|
&(currentAtomType->w0), |
| 593 |
|
&(currentAtomType->v0), |
| 594 |
+ |
&(currentAtomType->dipole), |
| 595 |
|
&isError ); |
| 596 |
|
if( isError ){ |
| 597 |
|
sprintf( painCave.errMsg, |
| 713 |
|
MPIcheckPoint(); |
| 714 |
|
#endif // is_mpi |
| 715 |
|
|
| 716 |
+ |
initFortran(); |
| 717 |
+ |
entry_plug->refreshSim(); |
| 718 |
+ |
|
| 719 |
|
} |
| 720 |
|
|
| 721 |
|
void TraPPE_ExFF::initializeBonds( bond_pair* the_bonds ){ |
| 1009 |
|
// if things go well, last will be set to 0 |
| 1010 |
|
|
| 1011 |
|
QuadraticBend* qBend; |
| 1012 |
+ |
GhostBend* gBend; |
| 1013 |
|
SRI **the_sris; |
| 1014 |
|
Atom** the_atoms; |
| 1015 |
|
int nBends; |
| 1110 |
|
|
| 1111 |
|
atomA = the_atoms[a]->getType(); |
| 1112 |
|
atomB = the_atoms[b]->getType(); |
| 1113 |
< |
atomC = the_atoms[c]->getType(); |
| 1113 |
> |
|
| 1114 |
> |
if( the_bends[i].isGhost ) atomC = "GHOST"; |
| 1115 |
> |
else atomC = the_atoms[c]->getType(); |
| 1116 |
> |
|
| 1117 |
|
currentBendType = headBendType->find( atomA, atomB, atomC ); |
| 1118 |
|
if( currentBendType == NULL ){ |
| 1119 |
|
sprintf( painCave.errMsg, "BendType error, %s - %s - %s not found" |
| 1126 |
|
if( !strcmp( currentBendType->type, "quadratic" ) ){ |
| 1127 |
|
|
| 1128 |
|
index = i + entry_plug->n_bonds; |
| 1129 |
< |
qBend = new QuadraticBend( *the_atoms[a], |
| 1130 |
< |
*the_atoms[b], |
| 1131 |
< |
*the_atoms[c] ); |
| 1132 |
< |
qBend->setConstants( currentBendType->k1, |
| 1133 |
< |
currentBendType->k2, |
| 1134 |
< |
currentBendType->k3, |
| 1135 |
< |
currentBendType->t0 ); |
| 1136 |
< |
the_sris[index] = qBend; |
| 1129 |
> |
|
| 1130 |
> |
if( the_bends[i].isGhost){ |
| 1131 |
> |
|
| 1132 |
> |
if( the_bends[i].ghost == b ){ |
| 1133 |
> |
// do nothing |
| 1134 |
> |
} |
| 1135 |
> |
else if( the_bends[i].ghost == a ){ |
| 1136 |
> |
c = a; |
| 1137 |
> |
a = b; |
| 1138 |
> |
b = a; |
| 1139 |
> |
} |
| 1140 |
> |
else{ |
| 1141 |
> |
sprintf( painCave.errMsg, |
| 1142 |
> |
"BendType error, %s - %s - %s,\n" |
| 1143 |
> |
" --> central atom is not " |
| 1144 |
> |
"correctly identified with the " |
| 1145 |
> |
"\"ghostVectorSource = \" tag.\n", |
| 1146 |
> |
atomA, atomB, atomC ); |
| 1147 |
> |
painCave.isFatal = 1; |
| 1148 |
> |
simError(); |
| 1149 |
> |
} |
| 1150 |
> |
|
| 1151 |
> |
gBend = new GhostBend( *the_atoms[a], |
| 1152 |
> |
*the_atoms[b] ); |
| 1153 |
> |
gBend->setConstants( currentBendType->k1, |
| 1154 |
> |
currentBendType->k2, |
| 1155 |
> |
currentBendType->k3, |
| 1156 |
> |
currentBendType->t0 ); |
| 1157 |
> |
the_sris[index] = gBend; |
| 1158 |
> |
} |
| 1159 |
> |
else{ |
| 1160 |
> |
qBend = new QuadraticBend( *the_atoms[a], |
| 1161 |
> |
*the_atoms[b], |
| 1162 |
> |
*the_atoms[c] ); |
| 1163 |
> |
qBend->setConstants( currentBendType->k1, |
| 1164 |
> |
currentBendType->k2, |
| 1165 |
> |
currentBendType->k3, |
| 1166 |
> |
currentBendType->t0 ); |
| 1167 |
> |
the_sris[index] = qBend; |
| 1168 |
> |
} |
| 1169 |
|
} |
| 1170 |
|
} |
| 1171 |
|
|
