ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/WATER.cpp
(Generate patch)

Comparing trunk/src/UseTheForce/WATER.cpp (file contents):
Revision 2 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
Revision 240 by tim, Fri Dec 10 18:41:05 2004 UTC

# Line 9 | Line 9 | using namespace std;
9   #include <mpi.h>
10   #endif //is_mpi
11  
12 < #include "ForceFields.hpp"
13 < #include "SRI.hpp"
14 < #include "simError.h"
12 > #include "UseTheForce/ForceFields.hpp"
13 > #include "primitives/SRI.hpp"
14 > #include "utils/simError.h"
15 > #include "types/AtomType.hpp"
16 > #include "types/DirectionalAtomType.hpp"
17 > #include "UseTheForce/DarkSide/lj_interface.h"
18 > #include "UseTheForce/DarkSide/charge_interface.h"
19 > #include "UseTheForce/DarkSide/dipole_interface.h"
20 > #include "UseTheForce/DarkSide/sticky_interface.h"
21  
16 – #include "fortranWrappers.hpp"
17 –
22   #ifdef IS_MPI
23 < #include "mpiForceField.h"
23 > #include "UseTheForce/mpiForceField.h"
24   #endif // is_mpi
25  
26  
# Line 243 | Line 247 | WATER::WATER(){
247  
248   WATER::WATER(){
249  
250 <  char fileName[200];
251 <  char* ffPath_env = "FORCE_PARAM_PATH";
248 <  char* ffPath;
249 <  char temp[200];
250 >  string fileName;
251 >  string tempString;
252  
253    headAtomType = NULL;
254    currentAtomType = NULL;
255    headDirectionalType = NULL;
256    currentDirectionalType = NULL;
257  
256 –  // do the funtion wrapping
257 –  wrapMeFF( this );
258 –
258   #ifdef IS_MPI
259    int i;
260    
# Line 309 | Line 308 | WATER::WATER(){
308    if( worldRank == 0 ){
309   #endif
310      
311 <    // generate the force file name
312 <    
313 <    strcpy( fileName, "WATER.frc" );
311 >    // generate the force file name  
312 >
313 >    fileName = "WATER.frc";
314 >
315      //    fprintf( stderr,"Trying to open %s\n", fileName );
316      
317      // attempt to open the file in the current directory first.
318      
319 <    frcFile = fopen( fileName, "r" );
319 >    frcFile = fopen( fileName.c_str(), "r" );
320      
321      if( frcFile == NULL ){
322        
323 <      // next see if the force path enviorment variable is set
323 >      tempString = ffPath + "/" + fileName;
324 >      fileName = tempString;
325        
326 <      ffPath = getenv( ffPath_env );
326 <      if( ffPath == NULL ) {
327 <        STR_DEFINE(ffPath, FRC_PATH );
328 <      }
326 >      frcFile = fopen( fileName.c_str(), "r" );
327        
330 –      
331 –      strcpy( temp, ffPath );
332 –      strcat( temp, "/" );
333 –      strcat( temp, fileName );
334 –      strcpy( fileName, temp );
335 –      
336 –      frcFile = fopen( fileName, "r" );
337 –      
328        if( frcFile == NULL ){
329          
330          sprintf( painCave.errMsg,
# Line 342 | Line 332 | WATER::WATER(){
332                   "\t%s\n"
333                   "\tHave you tried setting the FORCE_PARAM_PATH environment "
334                   "variable?\n",
335 <                 fileName );
335 >                 fileName.c_str() );
336          painCave.severity = OOPSE_ERROR;
337          painCave.isFatal = 1;
338          simError();
# Line 392 | Line 382 | void WATER::cleanMe( void ){
382   }
383  
384  
385 < void WATER::initForceField( int ljMixRule ){
385 > void WATER::initForceField(){
386    
387 <  initFortran( ljMixRule, entry_plug->useReactionField );
387 >  initFortran(entry_plug->useReactionField );
388   }
389  
390  
391   void WATER::readParams( void ){
392  
393 <  int identNum;
393 >  int identNum, isError;
394    int tempDirect0, tempDirect1;
395  
396    atomStruct atomInfo;
397    directionalStruct directionalInfo;
398    fpos_t *atomPos;
399 +  AtomType* at;
400  
401    atomInfo.last = 1;         // initialize last to have the last set.
402    directionalInfo.last = 1;  // if things go well, last will be set to 0
# Line 544 | Line 535 | void WATER::readParams( void ){
535    }
536  
537   #endif // is_mpi
547 –  
548 –  // call new A_types in fortran
538    
550 –  int isError;
551 –
539    // dummy variables
553 –  int isGB = 0;
554 –  int isEAM = 0;
555 –  int notDipole = 0;
556 –  int notSSD = 0;
557 –  double noDipMoment = 0.0;
540  
559 –
541    currentAtomType = headAtomType->next;
542    currentDirectionalType = headDirectionalType->next;
543  
544    while( currentAtomType != NULL ){
545 <
546 <    if( currentAtomType->isLJ ) entry_plug->useLJ = 1;
547 <    if( currentAtomType->isCharge ) entry_plug->useCharges = 1;
548 <    if( currentAtomType->isDirectional ){
549 <      // only directional atoms can have dipoles or be sticky
550 <      if ( currentDirectionalType->isDipole ) entry_plug->useDipoles = 1;
551 <      if ( currentDirectionalType->isSticky ) {
552 <        entry_plug->useSticky = 1;
572 <        set_sticky_params( &(currentDirectionalType->w0),
573 <                           &(currentDirectionalType->v0),
574 <                           &(currentDirectionalType->v0p),
575 <                           &(currentDirectionalType->rl),
576 <                           &(currentDirectionalType->ru),
577 <                           &(currentDirectionalType->rlp),
578 <                           &(currentDirectionalType->rup));
545 >    if( currentAtomType->name[0] != '\0' ){
546 >      if (currentAtomType->isDirectional)
547 >        at = new DirectionalAtomType();        
548 >      else
549 >        at = new AtomType();
550 >      
551 >      if (currentAtomType->isLJ) {
552 >        at->setLennardJones();
553        }
554 <      if( currentAtomType->name[0] != '\0' ){
555 <        isError = 0;
556 <        makeAtype( &(currentAtomType->ident),
557 <                   &(currentAtomType->isLJ),
558 <                   &(currentDirectionalType->isSticky),
559 <                   &(currentDirectionalType->isDipole),
560 <                   &isGB,
561 <                   &isEAM,
588 <                   &(currentAtomType->isCharge),
589 <                   &(currentAtomType->epslon),
590 <                   &(currentAtomType->sigma),
591 <                   &(currentAtomType->charge),
592 <                   &(currentDirectionalType->dipole),
593 <                   &isError );
594 <        if( isError ){
595 <          sprintf( painCave.errMsg,
596 <                   "Error initializing the \"%s\" atom type in fortran\n",
597 <                   currentAtomType->name );
598 <          painCave.isFatal = 1;
599 <          simError();
554 >
555 >      if (currentAtomType->isCharge) {
556 >        at->setCharge();
557 >      }
558 >
559 >      if (currentAtomType->isDirectional) {
560 >        if (currentDirectionalType->isDipole) {
561 >          ((DirectionalAtomType*)at)->setDipole();
562          }
563 +      
564 +        if (currentDirectionalType->isSticky) {
565 +          ((DirectionalAtomType*)at)->setSticky();
566 +        }
567        }
568 <      currentDirectionalType->next;
568 >      
569 >      at->setIdent(currentAtomType->ident);
570 >      at->setName(currentAtomType->name);    
571 >      at->complete();
572      }
573 +    currentAtomType = currentAtomType->next;
574 +  }
575  
576 <    else {
577 <      // use all dummy variables if this is not a directional atom
578 <      if( currentAtomType->name[0] != '\0' ){
579 <        isError = 0;
580 <        makeAtype( &(currentAtomType->ident),
581 <                   &(currentAtomType->isLJ),
582 <                   &notSSD,
583 <                   &notDipole,
584 <                   &isGB,
585 <                   &isEAM,
586 <                   &(currentAtomType->isCharge),
587 <                   &(currentAtomType->epslon),
588 <                   &(currentAtomType->sigma),
589 <                   &(currentAtomType->charge),
590 <                   &noDipMoment,
591 <                   &isError );
576 >  currentAtomType = headAtomType->next;
577 >  currentDirectionalType = headDirectionalType->next;
578 >
579 >  while( currentAtomType != NULL ){    
580 >
581 >    currentDirectionalType = headDirectionalType->find(currentAtomType->name);
582 >  
583 >    if( currentAtomType->isLJ ){
584 >      isError = 0;
585 >      newLJtype( &(currentAtomType->ident), &(currentAtomType->sigma),
586 >                 &(currentAtomType->epslon), &isError);
587 >      if( isError ){
588 >        sprintf( painCave.errMsg,
589 >                 "Error initializing the \"%s\" LJ type in fortran\n",
590 >                 currentAtomType->name );
591 >        painCave.isFatal = 1;
592 >        simError();
593 >      }
594 >    }
595 >
596 >    if (currentAtomType->isCharge) {
597 >      newChargeType(&(currentAtomType->ident), &(currentAtomType->charge),
598 >                    &isError);
599 >      if( isError ){
600 >        sprintf( painCave.errMsg,
601 >                 "Error initializing the \"%s\" charge type in fortran\n",
602 >                 currentAtomType->name );
603 >        painCave.isFatal = 1;
604 >        simError();
605 >      }
606 >    }
607 >
608 >    if (currentAtomType->isDirectional){
609 >      if (currentDirectionalType->isDipole) {
610 >        newDipoleType(&(currentAtomType->ident),
611 >                      &(currentDirectionalType->dipole),
612 >                      &isError);
613          if( isError ){
614            sprintf( painCave.errMsg,
615 <                   "Error initializing the \"%s\" atom type in fortran\n",
616 <                   currentAtomType->name );
615 >                   "Error initializing the \"%s\" dipole type in fortran\n",
616 >                   currentDirectionalType->name );
617            painCave.isFatal = 1;
618            simError();
619          }
620        }
621 +      
622 +      if(currentDirectionalType->isSticky) {        
623 +        makeStickyType( &(currentDirectionalType->w0),
624 +                        &(currentDirectionalType->v0),
625 +                        &(currentDirectionalType->v0p),
626 +                        &(currentDirectionalType->rl),
627 +                        &(currentDirectionalType->ru),
628 +                        &(currentDirectionalType->rlp),
629 +                        &(currentDirectionalType->rup));
630 +      }
631      }
632      currentAtomType = currentAtomType->next;
633    }
634 <
634 >  
635   #ifdef IS_MPI
636    sprintf( checkPointMsg,
637             "WATER atom and directional structures successfully"
638             "sent to fortran\n" );
639    MPIcheckPoint();
640   #endif // is_mpi
641 <
641 >  
642   }
643  
644  
# Line 658 | Line 660 | void WATER::initializeAtoms( int nAtoms, Atom** the_at
660        painCave.isFatal = 1;
661        simError();
662      }
661 –    the_atoms[i]->setMass( currentAtomType->mass );
662 –    the_atoms[i]->setIdent( currentAtomType->ident );
663 –
664 –    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
665 –
666 –    the_atoms[i]->setHasCharge(currentAtomType->isCharge);
663  
664      if( currentAtomType->isDirectional ){
665        currentDirectionalType =
# Line 676 | Line 672 | void WATER::initializeAtoms( int nAtoms, Atom** the_at
672          simError();
673        }
674  
675 +    the_atoms[i]->setMass( currentAtomType->mass );
676 +    the_atoms[i]->setIdent( currentAtomType->ident );
677 +
678 +    if (currentAtomType->isLJ) entry_plug->useLennardJones = 1;
679 +    if (currentAtomType->isCharge) entry_plug->useCharges = 1;
680 +    if (currentAtomType->isDirectional) {
681 +      if (currentDirectionalType->isDipole) entry_plug->useDipoles = 1;      
682 +      if (currentDirectionalType->isSticky) entry_plug->useSticky = 1;
683 +    }
684 +
685 +    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
686 +
687 +    the_atoms[i]->setHasCharge(currentAtomType->isCharge);
688 +
689 +
690 +
691        // zero out the moments of inertia matrix
692        for( j=0; j<3; j++ )
693          for( k=0; k<3; k++ )
# Line 1075 | Line 1087 | int WATER_NS::parseDirectional( char *lineBuffer, int
1087      return 1;
1088    }
1089    else return 0;
1078 – }
1079 – double WATER::getAtomTypeMass (char* atomType) {
1080 –
1081 –  currentAtomType = headAtomType->find( atomType );
1082 –  if( currentAtomType == NULL ){
1083 –    sprintf( painCave.errMsg,
1084 –            "AtomType error, %s not found in force file.\n",
1085 –             atomType );
1086 –    painCave.isFatal = 1;
1087 –    simError();
1088 –  }
1089 –
1090 –  return currentAtomType->mass;
1090   }

Diff Legend

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