ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/TraPPE_ExFF.cpp
(Generate patch)

Comparing trunk/OOPSE_old/src/mdtools/libmdCode/TraPPE_ExFF.cpp (file contents):
Revision 293 by mmeineke, Thu Mar 6 16:07:57 2003 UTC vs.
Revision 321 by mmeineke, Wed Mar 12 15:12:24 2003 UTC

# Line 9 | Line 9 | using namespace std;
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
# Line 84 | Line 86 | namespace TPE {  // restrict the access of the folowin
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  
# Line 108 | Line 103 | TraPPE_ExFF::TraPPE_ExFF(){
103    char errMsg[1000];
104  
105    // do the funtion wrapping
106 <  currentTPEwrap = this;
112 <  wrapMe();
106 >  wrapMeFF( this );
107  
108  
109   #ifdef IS_MPI
# Line 223 | Line 217 | TraPPE_ExFF::TraPPE_ExFF(){
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        
# Line 275 | Line 263 | TraPPE_ExFF::~TraPPE_ExFF(){
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 {
# Line 540 | Line 563 | void TraPPE_ExFF::initializeAtoms( void ){
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      
# Line 549 | Line 579 | void TraPPE_ExFF::initializeAtoms( void ){
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,
# Line 675 | Line 707 | void TraPPE_ExFF::initializeAtoms( void ){
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 ){
# Line 968 | Line 1003 | void TraPPE_ExFF::initializeBends( bend_set* the_bends
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;
# Line 1068 | Line 1104 | void TraPPE_ExFF::initializeBends( bend_set* the_bends
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"
# Line 1081 | Line 1120 | void TraPPE_ExFF::initializeBends( bend_set* the_bends
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  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines