mdtools/libmdCode/LJ_FF.cpp

#include <cstdlib>
#include <cstdio>
#include <cstring>

#include <iostream>
using namespace std;

#ifdef IS_MPI
#include <mpi.h>
#include <mpi++.h>
#endif //is_mpi

#include "ForceFields.hpp"
#include "SRI.hpp"
#include "simError.h"

#ifdef IS_MPI
#include "mpiForceField.h"
#endif // is_mpi


namespace {

  // Declare the structures that will be passed by the parser and  MPI
  
  typedef struct{
    char name[15];
    double mass;
    double epslon;
    double sigma;
    int ident;
    int last;      //  0  -> default
                   //  1  -> in MPI: tells nodes to stop listening
  } atomStruct;

  int parseAtom( char *lineBuffer, int lineNum, atomStruct &info );
  
#ifdef IS_MPI
  
  MPI_Datatype mpiAtomStructType;
  
#endif
}

// declaration of functions needed to wrap the fortran module

extern "C" {
  void forcefactory_( char* forceName,
                      int* status,
                      void (*wrapFunction)( void (*p1)( int* ident, 
                                                        double* mass, 
                                                        double* epslon, 
                                                        double* sigma, 
                                                        int* status ),
                                            void (*p2)( int *nLocal,
                                                        int *identArray,
                                                        int *isError ),
                                            void (*p3)( double* positionArray,
                                                        double* forceArray,
                                                        double* potentialEnergy,
                                                        double* tau,
                                                        short int* doPotentialCalc,
                                                        int* isError)),
                      int forceNameLength );
}


void LJfunctionWrapper( void (*p1)( int* ident, double* mass, double* epslon, 
                                   double* sigma, int* status ),
                        void (*p2)( int *nLocal, int *identArray, int *isError ),
                        void (*p3)( double* positionArray,double* forceArray,
                                    double* potentialEnergy, double* tau,
                                    short int* doPotentialCalc, int* isError ) );

void (*newLJtype)( int* ident, double* mass, double* epslon, double* sigma,
                   int* status );

void (*initLJfortran) ( int *nLocal, int *identArray, int *isError );

LJ_FF* currentLJwrap;


//****************************************************************
// begins the actual forcefield stuff.  
//****************************************************************


LJ_FF::LJ_FF(){

  char fileName[200];
  char* ffPath_env = "FORCE_PARAM_PATH";
  char* ffPath;
  char temp[200];
  char errMsg[1000];

  // do the funtion wrapping
  currentLJwrap = this;
  wrapMe();

#ifdef IS_MPI
  int i;
  
   // **********************************************************************
  // Init the atomStruct mpi type

  atomStruct atomProto; // mpiPrototype
  int atomBC[3] = {15,3,2};  // block counts
  MPI_Aint atomDspls[3];           // displacements
  MPI_Datatype atomMbrTypes[3];    // member mpi types

  MPI_Address(&atomProto.name, &atomDspls[0]);
  MPI_Address(&atomProto.mass, &atomDspls[1]);
  MPI_Address(&atomProto.ident, &atomDspls[2]);
  
  atomMbrTypes[0] = MPI_CHAR;
  atomMbrTypes[1] = MPI_DOUBLE;
  atomMbrTypes[2] = MPI_INT;
  
  for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
  
  MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
  MPI_Type_commit(&mpiAtomStructType);

  // ***********************************************************************
  
  if( worldRank == 0 ){
#endif
    
    // generate the force file name
    
    strcpy( fileName, "LJ_FF.frc" ); 
    //    fprintf( stderr,"Trying to open %s\n", fileName );
    
    // attempt to open the file in the current directory first.
    
    frcFile = fopen( fileName, "r" );
    
    if( frcFile == NULL ){
      
      // next see if the force path enviorment variable is set
      
      ffPath = getenv( ffPath_env );
      if( ffPath == NULL ) {
        sprintf( painCave.errMsg,
                 "Error opening the force field parameter file: %s\n"
                 "Have you tried setting the FORCE_PARAM_PATH environment "
                 "vairable?\n",
                 fileName );
        painCave.isFatal = 1;
        simError();
      }
      
      
      strcpy( temp, ffPath );
      strcat( temp, "/" );
      strcat( temp, fileName );
      strcpy( fileName, temp );
      
      frcFile = fopen( fileName, "r" );
      
      if( frcFile == NULL ){
        
        sprintf( painCave.errMsg,
                 "Error opening the force field parameter file: %s\n"
                 "Have you tried setting the FORCE_PARAM_PATH environment "
                 "vairable?\n",
                 fileName );
        painCave.isFatal = 1;
        simError();
      }
    }
    
#ifdef IS_MPI
  }
  
  sprintf( checkPointMsg, "LJ_FF file opened sucessfully." );
  MPIcheckPoint();
  
#endif // is_mpi
}


LJ_FF::~LJ_FF(){

#ifdef IS_MPI
  if( worldRank == 0 ){
#endif // is_mpi
    
    fclose( frcFile );
    
#ifdef IS_MPI
  }
#endif // is_mpi
}


void LJ_FF::wrapMe( void ){
  
  char* currentFF = "LJ";
  int isError = 0;
  
  forcefactory_( currentFF, &isError, LJfunctionWrapper, strlen(currentFF) );

  if( isError ){
    
    sprintf( painCave.errMsg,
             "LJ_FF error: an error was returned from fortran when the "
             "the functions were being wrapped.\n" );
    painCave.isFatal = 1;
    simError();
  }

#ifdef IS_MPI
  sprintf( checkPointMsg, "LJ_FF functions succesfully wrapped." );
  MPIcheckPoint();
#endif // is_mpi
}
  

void LJfunctionWrapper( void (*p1)( int* ident, double* mass, double* epslon, 
                                   double* sigma, int* status ),
                        void (*p2)( int*, int*, int* ),
                        void (*p3)( double*,double*,double*,double*,
                                    short int*, int* ) ){
  
  
  newLJtype = p1;
  initLJfortran = p2;
  currentLJwrap->setLJfortran( p3 );
}


void LJ_FF::initializeAtoms( void ){
  
  class LinkedType {
  public:
    LinkedType(){ 
      next = NULL;
      name[0] = '\0';
    }
    ~LinkedType(){ if( next != NULL ) delete next; }

    LinkedType* find(char* key){ 
      if( !strcmp(name, key) ) return this;
      if( next != NULL ) return next->find(key);
      return NULL;
    }
    

    void add( atomStruct &info ){
    
      // check for duplicates
      
      if( !strcmp( info.name, name ) ){
        sprintf( painCave.errMsg,
                 "Duplicate LJ atom type \"%s\" found in "
                 "the LJ_FF param file./n",
                 name );
        painCave.isFatal = 1;
        simError();
      }
      
      if( next != NULL ) next->add(info);
      else{
        next = new LinkedType();
        strcpy(next->name, info.name);
        next->mass     = info.mass;
        next->epslon   = info.epslon;
        next->sigma    = info.sigma;
        next->ident    = info.ident;
      }
    }
    

#ifdef IS_MPI
    
    void duplicate( atomStruct &info ){
      strcpy(info.name, name);
      info.mass     = mass;
      info.epslon   = epslon;
      info.sigma    = sigma;
      info.ident    = ident;
      info.last     = 0;
    }


#endif

    char name[15];
    double mass;
    double epslon;
    double sigma;
    int ident;
    LinkedType* next;
  };
  
  LinkedType* headAtomType; 
  LinkedType* currentAtomType;
  atomStruct info;
  info.last = 1; // initialize last to have the last set. 
                 // if things go well, last will be set to 0

  int i;
  int identNum;
  
  Atom** the_atoms;
  int nAtoms;
  the_atoms = entry_plug->atoms;
  nAtoms = entry_plug->n_atoms;
  
  
#ifdef IS_MPI
  if( worldRank == 0 ){
#endif
    
    // read in the atom types.

    headAtomType = new LinkedType;
    
    fastForward( "AtomTypes", "initializeAtoms" );

    // we are now at the AtomTypes section.
    
    eof_test = fgets( readLine, sizeof(readLine), frcFile );
    lineNum++;
    
    
    // read a line, and start parseing out the atom types 

    if( eof_test == NULL ){
      sprintf( painCave.errMsg, 
               "Error in reading Atoms from force file at line %d.\n",
               lineNum );
      painCave.isFatal = 1;
      simError();
    }
    
    identNum = 1;
    // stop reading at end of file, or at next section
    while( readLine[0] != '#' && eof_test != NULL ){

      // toss comment lines
      if( readLine[0] != '!' ){
        
        // the parser returns 0 if the line was blank
        if( parseAtomLJ( readLine, lineNum, info ) ){
          info.ident = identNum;
          headAtomType->add( info );;
          identNum++;
        }
      }
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;
    }

#ifdef IS_MPI
    
    // send out the linked list to all the other processes

    sprintf( checkPointMsg,
             "LJ_FF atom structures read successfully." );
    MPIcheckPoint();

    currentAtomType = headAtomType->next; //skip the first element who is a place holder.
    while( currentAtomType != NULL ){
      currentAtomType->duplicate( info );

 
      sendFrcStruct( &info, mpiAtomStructType );

      sprintf( checkPointMsg,
               "successfully sent lJ force type: \"%s\"\n",
               info.name );
      MPIcheckPoint();

      currentAtomType = currentAtomType->next;
    }
    info.last = 1;
    sendFrcStruct( &info, mpiAtomStructType );
    
  }

  else{
    
    // listen for node 0 to send out the force params
    
    MPIcheckPoint();

    headAtomType = new LinkedType;
    recieveFrcStruct( &info, mpiAtomStructType );
    
    while( !info.last ){


      headAtomType->add( info );
      
      MPIcheckPoint();

      recieveFrcStruct( &info, mpiAtomStructType );
    }
  }
#endif // is_mpi

  // call new A_types in fortran
  
  int isError;
  currentAtomType = headAtomType;
  while( currentAtomType != NULL ){
    
    if( currentAtomType->name[0] != '\0' ){
      isError = 0;
          newLJtype( &(currentAtomType->ident),
                 &(currentAtomType->mass),
                 &(currentAtomType->epslon),
                 &(currentAtomType->sigma),
                 &isError );
      if( isError ){
        sprintf( painCave.errMsg,
                 "Error initializing the \"%s\" atom type in fortran\n",
                 currentAtomType->name );
        painCave.isFatal = 1;
        simError();
      }
    }
    currentAtomType = currentAtomType->next;
  }
      
#ifdef IS_MPI
  sprintf( checkPointMsg,
           "LJ_FF atom structures successfully sent to fortran\n" );
  MPIcheckPoint();
#endif // is_mpi

  
  // initialize the atoms
  
  double bigSigma = 0.0;
  Atom* thisAtom;

  for( i=0; i<nAtoms; i++ ){
    
    currentAtomType = headAtomType->find( the_atoms[i]->getType() );
    if( currentAtomType == NULL ){
      sprintf( painCave.errMsg, 
               "AtomType error, %s not found in force file.\n",
               the_atoms[i]->getType() );
      painCave.isFatal = 1;
      simError();
    }
    
    the_atoms[i]->setMass( currentAtomType->mass );
    the_atoms[i]->setEpslon( currentAtomType->epslon );
    the_atoms[i]->setSigma( currentAtomType->sigma );
    the_atoms[i]->setIdent( currentAtomType->ident );
    the_atoms[i]->setLJ();

    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
  }

  
#ifdef IS_MPI
  double tempBig = bigSigma;
  MPI::COMM_WORLD.Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX );
#endif  //is_mpi

  //calc rCut and rList

  entry_plug->rCut = 2.5 * bigSigma;
  if(entry_plug->rCut > (entry_plug->box_x / 2.0)) entry_plug->rCut = entry_plug->box_x / 2.0;
  if(entry_plug->rCut > (entry_plug->box_y / 2.0)) entry_plug->rCut = entry_plug->box_y / 2.0;
  if(entry_plug->rCut > (entry_plug->box_z / 2.0)) entry_plug->rCut = entry_plug->box_z / 2.0;

  entry_plug->rList = entry_plug->rCut + 1.0;

  // clean up the memory
  
  delete headAtomType;

#ifdef IS_MPI
  sprintf( checkPointMsg, "LJ_FF atoms initialized succesfully" );
  MPIcheckPoint();
#endif // is_mpi

  initFortran();
  entry_plug->refreshSim();

}

void LJ_FF::initializeBonds( bond_pair* the_bonds ){
  
    if( entry_plug->n_bonds ){
      sprintf( painCave.errMsg,
               "LJ_FF does not support bonds.\n" );
      painCave.isFatal = 1;
      simError();
    }
#ifdef IS_MPI
  MPIcheckPoint();
#endif // is_mpi

}

void LJ_FF::initializeBends( bend_set* the_bends ){

    if( entry_plug->n_bends ){
      sprintf( painCave.errMsg,
               "LJ_FF does not support bends.\n" );
      painCave.isFatal = 1;
      simError();
    }
#ifdef IS_MPI
  MPIcheckPoint();
#endif // is_mpi

}

void LJ_FF::initializeTorsions( torsion_set* the_torsions ){

    if( entry_plug->n_torsions ){
      sprintf( painCave.errMsg,
               "LJ_FF does not support torsions.\n" );
      painCave.isFatal = 1;
      simError();
    }
#ifdef IS_MPI
  MPIcheckPoint();
#endif // is_mpi

}

void LJ_FF::doForces( int calcPot ){

  int i, isError;
  double* frc;
  double* pos;
  double* tau;
  short int passedCalcPot = (short int)calcPot;

  // forces are zeroed here, before any are acumulated.
  // NOTE: do not rezero the forces in Fortran.

  for(i=0; i<entry_plug->n_atoms; i++){
    entry_plug->atoms[i]->zeroForces();
  }

  frc = Atom::getFrcArray();
  pos = Atom::getPosArray();
  tau = entry_plug->tau;

  isError = 0;
  doLJfortran( pos, frc, &(entry_plug->lrPot), tau, &passedCalcPot, &isError );


  if( isError ){
    sprintf( painCave.errMsg,
             "Error returned from the fortran force calculation.\n" );
    painCave.isFatal = 1;
    simError();
  }

#ifdef IS_MPI
  sprintf( checkPointMsg,
           "returned from the force calculation.\n" );
  MPIcheckPoint();
#endif // is_mpi

}
  
void LJ_FF::initFortran( void ){
  
  int nLocal = entry_plug->n_atoms;
  int *ident;
  int isError;
  int i;

  ident = new int[nLocal];

  for(i=0; i<nLocal; i++){
    ident[i] = entry_plug->atoms[i]->getIdent();
  }

  isError = 0;
  initLJfortran( &nLocal, ident, &isError );
  
  if(isError){
    sprintf( painCave.errMsg,
             "LJ_FF error: There was an error initializing the component list in fortran.\n" );
    painCave.isFatal = 1;
    simError();
  }

  
#ifdef IS_MPI
  sprintf( checkPointMsg, "LJ_FF successfully initialized the fortran component list.\n" );
  MPIcheckPoint();
#endif // is_mpi
  
  delete[] ident;

}

  
void LJ_FF::fastForward( char* stopText, char* searchOwner ){

  int foundText = 0;
  char* the_token;

  rewind( frcFile );
  lineNum = 0;

  eof_test = fgets( readLine, sizeof(readLine), frcFile );
  lineNum++;
  if( eof_test == NULL ){
    sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
             " file is empty.\n",
             searchOwner );
    painCave.isFatal = 1;
    simError();
  }
  
  
  while( !foundText ){
    while( eof_test != NULL && readLine[0] != '#' ){
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;
    }
    if( eof_test == NULL ){
      sprintf( painCave.errMsg,
               "Error fast forwarding force file for %s at "
               "line %d: file ended unexpectedly.\n",
               searchOwner,
               lineNum );
      painCave.isFatal = 1;
      simError();
    }
    
    the_token = strtok( readLine, " ,;\t#\n" );
    foundText = !strcmp( stopText, the_token );
    
    if( !foundText ){
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;
      
      if( eof_test == NULL ){
        sprintf( painCave.errMsg, 
                 "Error fast forwarding force file for %s at "
                 "line %d: file ended unexpectedly.\n",
                 searchOwner,
                 lineNum );
        painCave.isFatal = 1;
        simError();
      } 
    }
  }  
}


int parseAtom( char *lineBuffer, int lineNum,  atomStruct &info ){

  char* the_token;
  
  the_token = strtok( lineBuffer, " \n\t,;" );
  if( the_token != NULL ){
    
    strcpy( info.name, the_token );
          
    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
      sprintf( painCave.errMsg, 
               "Error parseing AtomTypes: line %d\n", lineNum );
      painCave.isFatal = 1;
      simError();
    }
    
    info.mass = atof( the_token );
    
    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
      sprintf( painCave.errMsg, 
               "Error parseing AtomTypes: line %d\n", lineNum );
      painCave.isFatal = 1;
      simError();
    }
        
    info.epslon = atof( the_token );
          
    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
      sprintf( painCave.errMsg, 
               "Error parseing AtomTypes: line %d\n", lineNum );
      painCave.isFatal = 1;
      simError();
    }
        
    info.sigma = atof( the_token );
    
    return 1;
  }
  else return 0;
}
Revision:	291
Committed:	Wed Mar 5 20:35:54 2003 UTC (21 years, 6 months ago) by mmeineke
File size:	15631 byte(s)
Log Message:	overhauled TraPPE_Ex forcfield in C++
#	User	Rev	Content
1	mmeineke	270	#include <cstdlib>
2			#include <cstdio>
3			#include <cstring>
4
5			#include <iostream>
6			using namespace std;
7
8	mmeineke	275	#ifdef IS_MPI
9			#include <mpi.h>
10			#include <mpi++.h>
11			#endif //is_mpi
12
13	mmeineke	270	#include "ForceFields.hpp"
14			#include "SRI.hpp"
15			#include "simError.h"
16
17	mmeineke	291	#ifdef IS_MPI
18			#include "mpiForceField.h"
19			#endif // is_mpi
20	mmeineke	275
21
22	mmeineke	270
23	mmeineke	291	namespace {
24	mmeineke	270
25	mmeineke	291	// Declare the structures that will be passed by the parser and MPI
26
27			typedef struct{
28			char name[15];
29			double mass;
30			double epslon;
31			double sigma;
32			int ident;
33			int last; // 0 -> default
34			// 1 -> in MPI: tells nodes to stop listening
35			} atomStruct;
36	mmeineke	270
37	mmeineke	291	int parseAtom( char *lineBuffer, int lineNum, atomStruct &info );
38
39	mmeineke	270	#ifdef IS_MPI
40	mmeineke	291
41			MPI_Datatype mpiAtomStructType;
42
43	mmeineke	270	#endif
44	mmeineke	291	}
45	mmeineke	270
46			// declaration of functions needed to wrap the fortran module
47
48			extern "C" {
49			void forcefactory_( char* forceName,
50			int* status,
51			void (wrapFunction)( void (p1)( int* ident,
52			double* mass,
53			double* epslon,
54			double* sigma,
55			int* status ),
56			void (p2)( int nLocal,
57			int *identArray,
58			int *isError ),
59			void (p3)( double positionArray,
60			double* forceArray,
61			double* potentialEnergy,
62	chuckv	290	double* tau,
63			short int* doPotentialCalc,
64			int* isError)),
65	mmeineke	270	int forceNameLength );
66			}
67
68
69			void LJfunctionWrapper( void (p1)( int ident, double* mass, double* epslon,
70			double* sigma, int* status ),
71			void (p2)( int nLocal, int identArray, int isError ),
72			void (p3)( double positionArray,double* forceArray,
73	chuckv	290	double* potentialEnergy, double* tau,
74			short int* doPotentialCalc, int* isError ) );
75	mmeineke	270
76			void (newLJtype)( int ident, double* mass, double* epslon, double* sigma,
77			int* status );
78
79			void (initLJfortran) ( int nLocal, int identArray, int isError );
80
81			LJ_FF* currentLJwrap;
82
83
84			//****************************************************************
85			// begins the actual forcefield stuff.
86			//****************************************************************
87
88
89			LJ_FF::LJ_FF(){
90
91			char fileName[200];
92			char* ffPath_env = "FORCE_PARAM_PATH";
93			char* ffPath;
94			char temp[200];
95			char errMsg[1000];
96
97			// do the funtion wrapping
98			currentLJwrap = this;
99			wrapMe();
100
101			#ifdef IS_MPI
102			int i;
103
104			// **********************************************************************
105			// Init the atomStruct mpi type
106
107			atomStruct atomProto; // mpiPrototype
108			int atomBC[3] = {15,3,2}; // block counts
109			MPI_Aint atomDspls[3]; // displacements
110			MPI_Datatype atomMbrTypes[3]; // member mpi types
111
112			MPI_Address(&atomProto.name, &atomDspls[0]);
113			MPI_Address(&atomProto.mass, &atomDspls[1]);
114			MPI_Address(&atomProto.ident, &atomDspls[2]);
115
116			atomMbrTypes[0] = MPI_CHAR;
117			atomMbrTypes[1] = MPI_DOUBLE;
118			atomMbrTypes[2] = MPI_INT;
119
120			for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
121
122			MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
123			MPI_Type_commit(&mpiAtomStructType);
124
125			// ***********************************************************************
126
127			if( worldRank == 0 ){
128			#endif
129
130			// generate the force file name
131
132			strcpy( fileName, "LJ_FF.frc" );
133			// fprintf( stderr,"Trying to open %s\n", fileName );
134
135			// attempt to open the file in the current directory first.
136
137			frcFile = fopen( fileName, "r" );
138
139			if( frcFile == NULL ){
140
141			// next see if the force path enviorment variable is set
142
143			ffPath = getenv( ffPath_env );
144			if( ffPath == NULL ) {
145			sprintf( painCave.errMsg,
146			"Error opening the force field parameter file: %s\n"
147			"Have you tried setting the FORCE_PARAM_PATH environment "
148			"vairable?\n",
149			fileName );
150			painCave.isFatal = 1;
151			simError();
152			}
153
154
155			strcpy( temp, ffPath );
156			strcat( temp, "/" );
157			strcat( temp, fileName );
158			strcpy( fileName, temp );
159
160			frcFile = fopen( fileName, "r" );
161
162			if( frcFile == NULL ){
163
164			sprintf( painCave.errMsg,
165			"Error opening the force field parameter file: %s\n"
166			"Have you tried setting the FORCE_PARAM_PATH environment "
167			"vairable?\n",
168			fileName );
169			painCave.isFatal = 1;
170			simError();
171			}
172			}
173
174			#ifdef IS_MPI
175			}
176
177			sprintf( checkPointMsg, "LJ_FF file opened sucessfully." );
178			MPIcheckPoint();
179
180			#endif // is_mpi
181			}
182
183
184			LJ_FF::~LJ_FF(){
185
186			#ifdef IS_MPI
187			if( worldRank == 0 ){
188			#endif // is_mpi
189
190			fclose( frcFile );
191
192			#ifdef IS_MPI
193			}
194			#endif // is_mpi
195			}
196
197
198			void LJ_FF::wrapMe( void ){
199
200			char* currentFF = "LJ";
201			int isError = 0;
202
203			forcefactory_( currentFF, &isError, LJfunctionWrapper, strlen(currentFF) );
204
205			if( isError ){
206
207			sprintf( painCave.errMsg,
208			"LJ_FF error: an error was returned from fortran when the "
209			"the functions were being wrapped.\n" );
210			painCave.isFatal = 1;
211			simError();
212			}
213
214			#ifdef IS_MPI
215			sprintf( checkPointMsg, "LJ_FF functions succesfully wrapped." );
216			MPIcheckPoint();
217			#endif // is_mpi
218			}
219
220
221			void LJfunctionWrapper( void (p1)( int ident, double* mass, double* epslon,
222			double* sigma, int* status ),
223			void (p2)( int, int, int ),
224	mmeineke	291	void (p3)( double,double,double,double*,
225			short int, int ) ){
226	mmeineke	270
227
228			newLJtype = p1;
229			initLJfortran = p2;
230			currentLJwrap->setLJfortran( p3 );
231			}
232
233
234
235			void LJ_FF::initializeAtoms( void ){
236
237			class LinkedType {
238			public:
239			LinkedType(){
240			next = NULL;
241			name[0] = '\0';
242			}
243			~LinkedType(){ if( next != NULL ) delete next; }
244
245			LinkedType* find(char* key){
246			if( !strcmp(name, key) ) return this;
247			if( next != NULL ) return next->find(key);
248			return NULL;
249			}
250
251
252			void add( atomStruct &info ){
253
254			// check for duplicates
255
256			if( !strcmp( info.name, name ) ){
257			sprintf( painCave.errMsg,
258			"Duplicate LJ atom type \"%s\" found in "
259			"the LJ_FF param file./n",
260			name );
261			painCave.isFatal = 1;
262			simError();
263			}
264
265			if( next != NULL ) next->add(info);
266			else{
267			next = new LinkedType();
268			strcpy(next->name, info.name);
269			next->mass = info.mass;
270			next->epslon = info.epslon;
271			next->sigma = info.sigma;
272			next->ident = info.ident;
273			}
274			}
275
276
277			#ifdef IS_MPI
278
279			void duplicate( atomStruct &info ){
280			strcpy(info.name, name);
281			info.mass = mass;
282			info.epslon = epslon;
283			info.sigma = sigma;
284			info.ident = ident;
285			info.last = 0;
286			}
287
288
289			#endif
290
291			char name[15];
292			double mass;
293			double epslon;
294			double sigma;
295			int ident;
296			LinkedType* next;
297			};
298
299			LinkedType* headAtomType;
300			LinkedType* currentAtomType;
301			atomStruct info;
302			info.last = 1; // initialize last to have the last set.
303			// if things go well, last will be set to 0
304
305			int i;
306			int identNum;
307
308			Atom** the_atoms;
309			int nAtoms;
310			the_atoms = entry_plug->atoms;
311			nAtoms = entry_plug->n_atoms;
312
313
314			#ifdef IS_MPI
315			if( worldRank == 0 ){
316			#endif
317
318			// read in the atom types.
319
320			headAtomType = new LinkedType;
321
322			fastForward( "AtomTypes", "initializeAtoms" );
323
324			// we are now at the AtomTypes section.
325
326			eof_test = fgets( readLine, sizeof(readLine), frcFile );
327			lineNum++;
328
329
330			// read a line, and start parseing out the atom types
331
332			if( eof_test == NULL ){
333			sprintf( painCave.errMsg,
334			"Error in reading Atoms from force file at line %d.\n",
335			lineNum );
336			painCave.isFatal = 1;
337			simError();
338			}
339
340			identNum = 1;
341			// stop reading at end of file, or at next section
342			while( readLine[0] != '#' && eof_test != NULL ){
343
344			// toss comment lines
345			if( readLine[0] != '!' ){
346
347			// the parser returns 0 if the line was blank
348			if( parseAtomLJ( readLine, lineNum, info ) ){
349			info.ident = identNum;
350			headAtomType->add( info );;
351			identNum++;
352			}
353			}
354			eof_test = fgets( readLine, sizeof(readLine), frcFile );
355			lineNum++;
356			}
357
358			#ifdef IS_MPI
359
360			// send out the linked list to all the other processes
361
362			sprintf( checkPointMsg,
363			"LJ_FF atom structures read successfully." );
364			MPIcheckPoint();
365
366			currentAtomType = headAtomType->next; //skip the first element who is a place holder.
367			while( currentAtomType != NULL ){
368			currentAtomType->duplicate( info );
369
370
371
372			sendFrcStruct( &info, mpiAtomStructType );
373
374			sprintf( checkPointMsg,
375			"successfully sent lJ force type: \"%s\"\n",
376			info.name );
377			MPIcheckPoint();
378
379			currentAtomType = currentAtomType->next;
380			}
381			info.last = 1;
382			sendFrcStruct( &info, mpiAtomStructType );
383
384			}
385
386			else{
387
388			// listen for node 0 to send out the force params
389
390			MPIcheckPoint();
391
392			headAtomType = new LinkedType;
393			recieveFrcStruct( &info, mpiAtomStructType );
394
395			while( !info.last ){
396
397
398
399			headAtomType->add( info );
400
401			MPIcheckPoint();
402
403			recieveFrcStruct( &info, mpiAtomStructType );
404			}
405			}
406			#endif // is_mpi
407
408			// call new A_types in fortran
409
410			int isError;
411			currentAtomType = headAtomType;
412			while( currentAtomType != NULL ){
413
414			if( currentAtomType->name[0] != '\0' ){
415			isError = 0;
416			newLJtype( &(currentAtomType->ident),
417			&(currentAtomType->mass),
418			&(currentAtomType->epslon),
419			&(currentAtomType->sigma),
420			&isError );
421			if( isError ){
422			sprintf( painCave.errMsg,
423			"Error initializing the \"%s\" atom type in fortran\n",
424			currentAtomType->name );
425			painCave.isFatal = 1;
426			simError();
427			}
428			}
429			currentAtomType = currentAtomType->next;
430			}
431
432			#ifdef IS_MPI
433			sprintf( checkPointMsg,
434			"LJ_FF atom structures successfully sent to fortran\n" );
435			MPIcheckPoint();
436			#endif // is_mpi
437
438
439
440			// initialize the atoms
441
442			double bigSigma = 0.0;
443			Atom* thisAtom;
444
445			for( i=0; i<nAtoms; i++ ){
446
447			currentAtomType = headAtomType->find( the_atoms[i]->getType() );
448			if( currentAtomType == NULL ){
449			sprintf( painCave.errMsg,
450			"AtomType error, %s not found in force file.\n",
451			the_atoms[i]->getType() );
452			painCave.isFatal = 1;
453			simError();
454			}
455
456			the_atoms[i]->setMass( currentAtomType->mass );
457			the_atoms[i]->setEpslon( currentAtomType->epslon );
458			the_atoms[i]->setSigma( currentAtomType->sigma );
459			the_atoms[i]->setIdent( currentAtomType->ident );
460			the_atoms[i]->setLJ();
461
462			if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
463			}
464
465
466			#ifdef IS_MPI
467			double tempBig = bigSigma;
468			MPI::COMM_WORLD.Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX );
469			#endif //is_mpi
470
471			//calc rCut and rList
472
473			entry_plug->rCut = 2.5 * bigSigma;
474			if(entry_plug->rCut > (entry_plug->box_x / 2.0)) entry_plug->rCut = entry_plug->box_x / 2.0;
475			if(entry_plug->rCut > (entry_plug->box_y / 2.0)) entry_plug->rCut = entry_plug->box_y / 2.0;
476			if(entry_plug->rCut > (entry_plug->box_z / 2.0)) entry_plug->rCut = entry_plug->box_z / 2.0;
477
478			entry_plug->rList = entry_plug->rCut + 1.0;
479
480			// clean up the memory
481
482			delete headAtomType;
483
484			#ifdef IS_MPI
485			sprintf( checkPointMsg, "LJ_FF atoms initialized succesfully" );
486			MPIcheckPoint();
487			#endif // is_mpi
488
489			initFortran();
490			entry_plug->refreshSim();
491
492			}
493
494			void LJ_FF::initializeBonds( bond_pair* the_bonds ){
495
496			if( entry_plug->n_bonds ){
497			sprintf( painCave.errMsg,
498			"LJ_FF does not support bonds.\n" );
499			painCave.isFatal = 1;
500			simError();
501			}
502			#ifdef IS_MPI
503			MPIcheckPoint();
504			#endif // is_mpi
505
506			}
507
508			void LJ_FF::initializeBends( bend_set* the_bends ){
509
510			if( entry_plug->n_bends ){
511			sprintf( painCave.errMsg,
512			"LJ_FF does not support bends.\n" );
513			painCave.isFatal = 1;
514			simError();
515			}
516			#ifdef IS_MPI
517			MPIcheckPoint();
518			#endif // is_mpi
519
520			}
521
522			void LJ_FF::initializeTorsions( torsion_set* the_torsions ){
523
524			if( entry_plug->n_torsions ){
525			sprintf( painCave.errMsg,
526			"LJ_FF does not support torsions.\n" );
527			painCave.isFatal = 1;
528			simError();
529			}
530			#ifdef IS_MPI
531			MPIcheckPoint();
532			#endif // is_mpi
533
534			}
535
536	mmeineke	291	void LJ_FF::doForces( int calcPot ){
537	mmeineke	270
538	mmeineke	291	int i, isError;
539			double* frc;
540			double* pos;
541			double* tau;
542			short int passedCalcPot = (short int)calcPot;
543
544			// forces are zeroed here, before any are acumulated.
545			// NOTE: do not rezero the forces in Fortran.
546
547			for(i=0; i<entry_plug->n_atoms; i++){
548			entry_plug->atoms[i]->zeroForces();
549			}
550
551			frc = Atom::getFrcArray();
552			pos = Atom::getPosArray();
553			tau = entry_plug->tau;
554
555			isError = 0;
556			doLJfortran( pos, frc, &(entry_plug->lrPot), tau, &passedCalcPot, &isError );
557
558
559			if( isError ){
560			sprintf( painCave.errMsg,
561			"Error returned from the fortran force calculation.\n" );
562			painCave.isFatal = 1;
563			simError();
564			}
565
566			#ifdef IS_MPI
567			sprintf( checkPointMsg,
568			"returned from the force calculation.\n" );
569			MPIcheckPoint();
570			#endif // is_mpi
571
572			}
573
574			void LJ_FF::initFortran( void ){
575
576			int nLocal = entry_plug->n_atoms;
577			int *ident;
578			int isError;
579			int i;
580
581			ident = new int[nLocal];
582
583			for(i=0; i<nLocal; i++){
584			ident[i] = entry_plug->atoms[i]->getIdent();
585			}
586
587			isError = 0;
588			initLJfortran( &nLocal, ident, &isError );
589
590			if(isError){
591			sprintf( painCave.errMsg,
592			"LJ_FF error: There was an error initializing the component list in fortran.\n" );
593			painCave.isFatal = 1;
594			simError();
595			}
596
597
598			#ifdef IS_MPI
599			sprintf( checkPointMsg, "LJ_FF successfully initialized the fortran component list.\n" );
600			MPIcheckPoint();
601			#endif // is_mpi
602
603			delete[] ident;
604
605			}
606
607
608	mmeineke	270	void LJ_FF::fastForward( char* stopText, char* searchOwner ){
609
610			int foundText = 0;
611			char* the_token;
612
613			rewind( frcFile );
614			lineNum = 0;
615
616			eof_test = fgets( readLine, sizeof(readLine), frcFile );
617			lineNum++;
618			if( eof_test == NULL ){
619			sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
620			" file is empty.\n",
621			searchOwner );
622			painCave.isFatal = 1;
623			simError();
624			}
625
626
627			while( !foundText ){
628			while( eof_test != NULL && readLine[0] != '#' ){
629			eof_test = fgets( readLine, sizeof(readLine), frcFile );
630			lineNum++;
631			}
632			if( eof_test == NULL ){
633			sprintf( painCave.errMsg,
634			"Error fast forwarding force file for %s at "
635			"line %d: file ended unexpectedly.\n",
636			searchOwner,
637			lineNum );
638			painCave.isFatal = 1;
639			simError();
640			}
641
642			the_token = strtok( readLine, " ,;\t#\n" );
643			foundText = !strcmp( stopText, the_token );
644
645			if( !foundText ){
646			eof_test = fgets( readLine, sizeof(readLine), frcFile );
647			lineNum++;
648
649			if( eof_test == NULL ){
650			sprintf( painCave.errMsg,
651			"Error fast forwarding force file for %s at "
652			"line %d: file ended unexpectedly.\n",
653			searchOwner,
654			lineNum );
655			painCave.isFatal = 1;
656			simError();
657			}
658			}
659			}
660			}
661
662
663
664	mmeineke	291	int parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
665	mmeineke	270
666			char* the_token;
667
668			the_token = strtok( lineBuffer, " \n\t,;" );
669			if( the_token != NULL ){
670
671			strcpy( info.name, the_token );
672
673			if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
674			sprintf( painCave.errMsg,
675			"Error parseing AtomTypes: line %d\n", lineNum );
676			painCave.isFatal = 1;
677			simError();
678			}
679
680			info.mass = atof( the_token );
681
682			if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
683			sprintf( painCave.errMsg,
684			"Error parseing AtomTypes: line %d\n", lineNum );
685			painCave.isFatal = 1;
686			simError();
687			}
688
689			info.epslon = atof( the_token );
690
691			if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
692			sprintf( painCave.errMsg,
693			"Error parseing AtomTypes: line %d\n", lineNum );
694			painCave.isFatal = 1;
695			simError();
696			}
697
698			info.sigma = atof( the_token );
699
700			return 1;
701			}
702			else return 0;
703			}