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 296 by mmeineke, Thu Mar 6 20:05:39 2003 UTC vs.
Revision 373 by mmeineke, Thu Mar 20 19:01:44 2003 UTC

# Line 12 | Line 12 | using namespace std;
12   #include <fortranWrappers.hpp>
13  
14   #ifdef IS_MPI
15 + #include <mpi++.h>
16   #include "mpiForceField.h"
17   #endif // is_mpi
18  
# Line 217 | Line 218 | TraPPE_ExFF::TraPPE_ExFF(){
218        
219        ffPath = getenv( ffPath_env );
220        if( ffPath == NULL ) {
221 <        sprintf( painCave.errMsg,
221 <                 "Error opening the force field parameter file: %s\n"
222 <                 "Have you tried setting the FORCE_PARAM_PATH environment "
223 <                 "vairable?\n",
224 <                 fileName );
225 <        painCave.isFatal = 1;
226 <        simError();
221 >        STR_DEFINE(ffPath, FRC_PATH );
222        }
223        
224        
# Line 268 | Line 263 | TraPPE_ExFF::~TraPPE_ExFF(){
263    }
264   #endif // is_mpi
265   }
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;
266  
267 <  // 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 ){
267 > void TraPPE_ExFF::initForceField( int ljMixRule ){
268    
269 <  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 <
269 >  initFortran( ljMixRule, entry_plug->useReactionField );
270   }
271  
272  
345 –
346 –
273   void TraPPE_ExFF::initializeAtoms( void ){
274    
275    class LinkedType {
# Line 436 | Line 362 | void TraPPE_ExFF::initializeAtoms( void ){
362    the_atoms = entry_plug->atoms;
363    nAtoms = entry_plug->n_atoms;
364    
439 –  
365    //////////////////////////////////////////////////
366    // a quick water fix
367  
# Line 574 | Line 499 | void TraPPE_ExFF::initializeAtoms( void ){
499    
500    int isGB = 0;
501    int isLJ = 1;
502 +  double GB_dummy = 0.0;
503    
504    
505    currentAtomType = headAtomType;
506    while( currentAtomType != NULL ){
507      
508 +    if(currentAtomType->isDipole) entry_plug->useReactionField = 1;
509 +    if(currentAtomType->isDipole) entry_plug->useDipole = 1;
510 +    if(currentAtomType->isSSD)    entry_plug->useSticky = 1;
511 +
512      if( currentAtomType->name[0] != '\0' ){
513        isError = 0;
514 <          newTPEtype( &(currentAtomType->ident),
515 <                      &(currentAtomType->mass),
516 <                      &(currentAtomType->epslon),
517 <                      &(currentAtomType->sigma),
518 <                      &isLJ,
519 <                      &(currentAtomType->isSSD),
520 <                      &(currentAtomType->isDipole),
521 <                      &isGB,
522 <                      &(currentAtomType->w0),
523 <                      &(currentAtomType->v0),
524 <                      &(currentAtomType->dipole),
525 <                      &isError );
514 >      makeAtype( &(currentAtomType->ident),
515 >                 &isLJ,
516 >                 &(currentAtomType->isSSD),
517 >                 &(currentAtomType->isDipole),
518 >                 &isGB,
519 >                 &(currentAtomType->epslon),
520 >                 &(currentAtomType->sigma),
521 >                 &(currentAtomType->dipole),
522 >                 &(currentAtomType->w0),
523 >                 &(currentAtomType->v0),
524 >                 &GB_dummy,
525 >                 &GB_dummy,
526 >                 &GB_dummy,
527 >                 &GB_dummy,
528 >                 &GB_dummy,
529 >                 &GB_dummy,
530 >                 &isError );
531        if( isError ){
532          sprintf( painCave.errMsg,
533                   "Error initializing the \"%s\" atom type in fortran\n",
# Line 703 | Line 638 | void TraPPE_ExFF::initializeAtoms( void ){
638  
639    entry_plug->rList = entry_plug->rCut + 1.0;
640  
641 +  entry_plug->useLJ = 1; // use Lennard Jones is on by default
642  
643    // clean up the memory
644    
# Line 713 | Line 649 | void TraPPE_ExFF::initializeAtoms( void ){
649    MPIcheckPoint();
650   #endif // is_mpi
651  
716 –  initFortran();
717 –  entry_plug->refreshSim();
718 –
652   }
653  
654   void TraPPE_ExFF::initializeBonds( bond_pair* the_bonds ){
# Line 1009 | Line 942 | void TraPPE_ExFF::initializeBends( bend_set* the_bends
942                   // if things go well, last will be set to 0
943  
944    QuadraticBend* qBend;
945 +  GhostBend* gBend;
946    SRI **the_sris;
947    Atom** the_atoms;
948    int nBends;
# Line 1109 | Line 1043 | void TraPPE_ExFF::initializeBends( bend_set* the_bends
1043  
1044      atomA = the_atoms[a]->getType();
1045      atomB = the_atoms[b]->getType();
1046 <    atomC = the_atoms[c]->getType();
1046 >
1047 >    if( the_bends[i].isGhost ) atomC = "GHOST";
1048 >    else atomC = the_atoms[c]->getType();
1049 >
1050      currentBendType = headBendType->find( atomA, atomB, atomC );
1051      if( currentBendType == NULL ){
1052        sprintf( painCave.errMsg, "BendType error, %s - %s - %s not found"
# Line 1122 | Line 1059 | void TraPPE_ExFF::initializeBends( bend_set* the_bends
1059      if( !strcmp( currentBendType->type, "quadratic" ) ){
1060        
1061        index = i + entry_plug->n_bonds;
1062 <      qBend = new QuadraticBend( *the_atoms[a],
1063 <                                 *the_atoms[b],
1064 <                                 *the_atoms[c] );
1065 <      qBend->setConstants( currentBendType->k1,
1066 <                           currentBendType->k2,
1067 <                           currentBendType->k3,
1068 <                           currentBendType->t0 );
1069 <      the_sris[index] = qBend;
1062 >      
1063 >      if( the_bends[i].isGhost){
1064 >        
1065 >        if( the_bends[i].ghost == b ){
1066 >          // do nothing
1067 >        }
1068 >        else if( the_bends[i].ghost == a ){
1069 >          c = a;
1070 >          a = b;
1071 >          b = a;
1072 >        }
1073 >        else{
1074 >          sprintf( painCave.errMsg,
1075 >                   "BendType error, %s - %s - %s,\n"
1076 >                   "  --> central atom is not "
1077 >                   "correctly identified with the "
1078 >                   "\"ghostVectorSource = \" tag.\n",
1079 >                   atomA, atomB, atomC );
1080 >          painCave.isFatal = 1;
1081 >          simError();
1082 >        }
1083 >        
1084 >        gBend = new GhostBend( *the_atoms[a],
1085 >                               *the_atoms[b] );                        
1086 >        gBend->setConstants( currentBendType->k1,
1087 >                             currentBendType->k2,
1088 >                             currentBendType->k3,
1089 >                             currentBendType->t0 );
1090 >        the_sris[index] = gBend;
1091 >      }
1092 >      else{
1093 >        qBend = new QuadraticBend( *the_atoms[a],
1094 >                                   *the_atoms[b],
1095 >                                   *the_atoms[c] );
1096 >        qBend->setConstants( currentBendType->k1,
1097 >                             currentBendType->k2,
1098 >                             currentBendType->k3,
1099 >                             currentBendType->t0 );
1100 >        the_sris[index] = qBend;
1101 >      }
1102      }
1103    }
1104  
# Line 1159 | Line 1128 | void TraPPE_ExFF::initializeTorsions( torsion_set* the
1128      ~LinkedType(){ if( next != NULL ) delete next; }
1129  
1130      LinkedType* find( char* key1, char* key2, char* key3, char* key4 ){
1131 +      
1132 +      std::cerr<< "looking for: " << key1 << " - " << key2 << " - "
1133 +                  << key3 << " - " << key4 << "\n";
1134 +
1135 +      std::cerr<< "I've got: " << nameA << " - " << nameB << " - "
1136 +                  << nameC << " - " << nameD << "\n";
1137 +
1138 +
1139 +
1140        if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) &&
1141            !strcmp( nameC, key3 ) && !strcmp( nameD, key4 ) ) return this;
1142  
# Line 1194 | Line 1172 | void TraPPE_ExFF::initializeTorsions( torsion_set* the
1172          strcpy(next->nameA, info.nameA);
1173          strcpy(next->nameB, info.nameB);
1174          strcpy(next->nameC, info.nameC);
1175 +        strcpy(next->nameD, info.nameD);
1176          strcpy(next->type,  info.type);
1177          next->k1 = info.k1;
1178          next->k2 = info.k2;
1179          next->k3 = info.k3;
1180          next->k4 = info.k4;
1181 +
1182 +        std::cerr << "added: " << info.nameA << " - " << info.nameB << " - "
1183 +                  << info.nameC << " - " << info.nameD << "\n";
1184        }
1185      }
1186  
# Line 1284 | Line 1266 | void TraPPE_ExFF::initializeTorsions( torsion_set* the
1266          // the parser returns 0 if the line was blank
1267          if( parseTorsion( readLine, lineNum, info ) ){
1268            headTorsionType->add( info );
1269 +
1270          }
1271        }
1272        eof_test = fgets( readLine, sizeof(readLine), frcFile );
# Line 1426 | Line 1409 | int parseAtom( char *lineBuffer, int lineNum, atomStru
1409   }
1410  
1411  
1412 < int parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
1412 > int TPE::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
1413  
1414    char* the_token;
1415    
# Line 1520 | Line 1503 | int parseBond( char *lineBuffer, int lineNum, bondStru
1503    else return 0;
1504   }
1505  
1506 < int parseBond( char *lineBuffer, int lineNum, bondStruct &info ){
1506 > int TPE::parseBond( char *lineBuffer, int lineNum, bondStruct &info ){
1507  
1508    char* the_token;
1509    
# Line 1572 | Line 1555 | int parseBend( char *lineBuffer, int lineNum, bendStru
1555   }
1556  
1557  
1558 < int parseBend( char *lineBuffer, int lineNum, bendStruct &info ){
1558 > int TPE::parseBend( char *lineBuffer, int lineNum, bendStruct &info ){
1559  
1560    char* the_token;
1561    
# Line 1660 | Line 1643 | int parseTorsion( char *lineBuffer, int lineNum, torsi
1643    else return 0;
1644   }
1645  
1646 < int parseTorsion( char *lineBuffer, int lineNum, torsionStruct &info ){
1646 > int TPE::parseTorsion( char *lineBuffer, int lineNum, torsionStruct &info ){
1647    
1648    char*  the_token;
1649  

Diff Legend

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