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 294 by mmeineke, Thu Mar 6 17:04:09 2003 UTC vs.
Revision 311 by mmeineke, Tue Mar 11 16:27:34 2003 UTC

# Line 264 | Line 264 | TraPPE_ExFF::~TraPPE_ExFF(){
264      
265      fclose( frcFile );
266      
267 + #ifdef IS_MPI
268 +  }
269 + #endif // is_mpi
270 + }
271 +
272 + void TraPPE_ExFF::doForces( int calcPot ){
273 +
274 +  int i, isError;
275 +  double* frc;
276 +  double* pos;
277 +  double* tau;
278 +  short int passedCalcPot = (short int)calcPot;
279 +
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 +             "Error returned from the fortran force calculation.\n" );
299 +    painCave.isFatal = 1;
300 +    simError();
301 +  }
302 +
303   #ifdef IS_MPI
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 +  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 + #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 {
# Line 495 | Line 569 | void TraPPE_ExFF::initializeAtoms( void ){
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      
# Line 504 | Line 585 | void TraPPE_ExFF::initializeAtoms( void ){
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,
# Line 629 | Line 712 | void TraPPE_ExFF::initializeAtoms( void ){
712    sprintf( checkPointMsg, "TraPPE_Ex atoms initialized succesfully" );
713    MPIcheckPoint();
714   #endif // is_mpi
715 +
716 +  initFortran();
717 +  entry_plug->refreshSim();
718  
719   }
720  
# Line 923 | Line 1009 | void TraPPE_ExFF::initializeBends( bend_set* the_bends
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;
# Line 1023 | Line 1110 | void TraPPE_ExFF::initializeBends( bend_set* the_bends
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"
# Line 1036 | Line 1126 | void TraPPE_ExFF::initializeBends( bend_set* the_bends
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  

Diff Legend

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