--- trunk/OOPSE/libmdtools/WATER.cpp 2004/01/22 17:34:20 976 +++ trunk/OOPSE/libmdtools/WATER.cpp 2004/04/12 20:32:20 1097 @@ -23,7 +23,7 @@ namespace WATER_NS{ namespace WATER_NS{ - // Declare the structures that will be passed by the parser and MPI + // Declare the structures that will be passed by the parser and MPI typedef struct{ char name[15]; @@ -41,6 +41,9 @@ namespace WATER_NS{ typedef struct{ char name[15]; + double Ixx; + double Iyy; + double Izz; double dipole; double w0; double v0; @@ -172,7 +175,10 @@ namespace WATER_NS{ strcpy(next->name, info.name); next->isDipole = info.isDipole; next->isSticky = info.isSticky; - next->dipole = info.dipole; + next->Ixx = info.Ixx; + next->Iyy = info.Iyy; + next->Izz = info.Izz; + next->dipole = info.dipole; next->w0 = info.w0; next->v0 = info.v0; next->v0p = info.v0p; @@ -189,6 +195,9 @@ namespace WATER_NS{ strcpy(info.name, name); info.isDipole = isDipole; info.isSticky = isSticky; + info.Ixx = Ixx; + info.Iyy = Iyy; + info.Izz = Izz; info.dipole = dipole; info.w0 = w0; info.v0 = v0; @@ -205,6 +214,9 @@ namespace WATER_NS{ char name[15]; int isDipole; int isSticky; + double Ixx; + double Iyy; + double Izz; double dipole; double w0; double v0; @@ -220,7 +232,7 @@ namespace WATER_NS{ LinkedAtomType* currentAtomType; LinkedDirectionalType* headDirectionalType; LinkedDirectionalType* currentDirectionalType; -} +} // namespace using namespace WATER_NS; @@ -257,7 +269,7 @@ WATER::WATER(){ MPI_Address(&atomProto.name, &atomDspls[0]); MPI_Address(&atomProto.mass, &atomDspls[1]); - MPI_Address(&atomProto.ident, &atomDspls[2]); + MPI_Address(&atomProto.isDirectional, &atomDspls[2]); atomMbrTypes[0] = MPI_CHAR; atomMbrTypes[1] = MPI_DOUBLE; @@ -274,13 +286,13 @@ WATER::WATER(){ // Init the directionalStruct mpi type directionalStruct directionalProto; // mpiPrototype - int directionalBC[3] = {15,8,3}; // block counts + int directionalBC[3] = {15,11,3}; // block counts MPI_Aint directionalDspls[3]; // displacements MPI_Datatype directionalMbrTypes[3]; // member mpi types MPI_Address(&directionalProto.name, &directionalDspls[0]); - MPI_Address(&directionalProto.mass, &directionalDspls[1]); - MPI_Address(&directionalProto.ident, &directionalDspls[2]); + MPI_Address(&directionalProto.Ixx, &directionalDspls[1]); + MPI_Address(&directionalProto.isDipole, &directionalDspls[2]); directionalMbrTypes[0] = MPI_CHAR; directionalMbrTypes[1] = MPI_DOUBLE; @@ -293,7 +305,7 @@ WATER::WATER(){ MPI_Type_commit(&mpiDirectionalStructType); // *********************************************************************** - + if( worldRank == 0 ){ #endif @@ -387,16 +399,18 @@ void WATER::readParams( void ){ void WATER::readParams( void ){ int identNum; + int tempDirect0, tempDirect1; atomStruct atomInfo; directionalStruct directionalInfo; + fpos_t *atomPos; atomInfo.last = 1; // initialize last to have the last set. directionalInfo.last = 1; // if things go well, last will be set to 0 atomPos = new fpos_t; bigSigma = 0.0; - + #ifdef IS_MPI if( worldRank == 0 ){ #endif @@ -404,9 +418,10 @@ void WATER::readParams( void ){ // read in the atom types. headAtomType = new LinkedAtomType; + headDirectionalType = new LinkedDirectionalType; fastForward( "AtomTypes", "initializeAtoms" ); - + // we are now at the AtomTypes section. eof_test = fgets( readLine, sizeof(readLine), frcFile ); @@ -414,7 +429,7 @@ void WATER::readParams( void ){ // read a line, and start parsing out the atom types - + if( eof_test == NULL ){ sprintf( painCave.errMsg, "Error in reading Atoms from force file at line %d.\n", @@ -442,6 +457,7 @@ void WATER::readParams( void ){ sectionSearch( "DirectionalTypes", atomInfo.name, "initializeDirectional" ); parseDirectional( readLine, lineNum, directionalInfo ); + headDirectionalType->add( directionalInfo ); // return to the AtomTypes section fsetpos( frcFile, atomPos ); @@ -454,87 +470,78 @@ void WATER::readParams( void ){ } #ifdef IS_MPI - + // send out the linked list to all the other processes sprintf( checkPointMsg, - "WATER atom structures read successfully." ); + "WATER atom and directional structures read successfully." ); MPIcheckPoint(); + currentAtomType = headAtomType->next; //skip the first element place holder + currentDirectionalType = headDirectionalType->next; // same w/ directional - currentAtomType = headAtomType->next; //skip the first element who is a place holder. while( currentAtomType != NULL ){ currentAtomType->duplicate( atomInfo ); - + sendFrcStruct( &atomInfo, mpiAtomStructType ); sprintf( checkPointMsg, "successfully sent WATER force type: \"%s\"\n", atomInfo.name ); - MPIcheckPoint(); + + if ( atomInfo.isDirectional ){ + // send out the directional linked list to all the other processes + + currentDirectionalType->duplicate( directionalInfo ); + sendFrcStruct( &directionalInfo, mpiDirectionalStructType ); + + sprintf( checkPointMsg, + "successfully sent WATER directional type: \"%s\"\n", + directionalInfo.name ); + } + MPIcheckPoint(); + tempDirect0 = atomInfo.isDirectional; currentAtomType = currentAtomType->next; + if( tempDirect0 ) + currentDirectionalType = currentDirectionalType->next; } + atomInfo.last = 1; sendFrcStruct( &atomInfo, mpiAtomStructType ); - - if ( atomInfo.isDirectional ){ - // send out the linked list to all the other processes - - sprintf( checkPointMsg, - "WATER directional structures read successfully." ); - MPIcheckPoint(); - - currentDirectionalType = headDirectionalType->next; - while( currentDirectionalType != NULL ){ - currentDirectionalType->duplicate( directionalInfo ); - sendFrcStruct( &directionalInfo, mpiDirectionalStructType ); - currentDirectionalType = currentDirectionalType->next; - } - directionalInfo.last = 1; + directionalInfo.last = 1; + if ( atomInfo.isDirectional ) sendFrcStruct( &directionalInfo, mpiDirectionalStructType ); - } } else{ - // listen for node 0 to send out the force params MPIcheckPoint(); - + headAtomType = new LinkedAtomType; + headDirectionalType = new LinkedDirectionalType; receiveFrcStruct( &atomInfo, mpiAtomStructType ); - + + if ( atomInfo.isDirectional ) + receiveFrcStruct( &directionalInfo, mpiDirectionalStructType ); + while( !atomInfo.last ){ - + headAtomType->add( atomInfo ); MPIcheckPoint(); - + receiveFrcStruct( &atomInfo, mpiAtomStructType ); - } - if ( atomInfo.isDirectional ) { - // listen for node 0 to send out the force params - - MPIcheckPoint(); - - headDirectionalType = new LinkedDirectionalType; - receiveFrcStruct( &directionalInfo, mpiDirectionalStructType ); - while( !directionalInfo.last ){ - + if( atomInfo.isDirectional ){ headDirectionalType->add( directionalInfo ); + receiveFrcStruct( &directionalInfo, mpiDirectionalStructType ); } - - sprintf( checkPointMsg, - "WATER directional structures broadcast successfully." ); - MPIcheckPoint(); } } #endif // is_mpi - - // call new A_types in fortran @@ -543,43 +550,79 @@ void WATER::readParams( void ){ // dummy variables int isGB = 0; int isEAM = 0; - - currentAtomType = headAtomType; + int notDipole = 0; + int notSSD = 0; + double noDipMoment = 0.0; + + + currentAtomType = headAtomType->next; + currentDirectionalType = headDirectionalType->next; + while( currentAtomType != NULL ){ if( currentAtomType->isLJ ) entry_plug->useLJ = 1; - if( currentAtomType->isCharge ) entry_plug->useCharge = 1; + if( currentAtomType->isCharge ) entry_plug->useCharges = 1; if( currentAtomType->isDirectional ){ - if ( currentAtomType->isDipole ){ - entry_plug->useDipoles = 1; - &(currentAtomType->dipole); - } - if ( currentAtomType->isSticky ) { + // only directional atoms can have dipoles or be sticky + if ( currentDirectionalType->isDipole ) entry_plug->useDipoles = 1; + if ( currentDirectionalType->isSticky ) { entry_plug->useSticky = 1; - set_sticky_params( &(currentAtomType->w0), &(currentAtomType->v0), - &(currentAtomType->v0p), - &(currentAtomType->rl), &(currentAtomType->ru), - &(currentAtomType->rlp), &(currentAtomType->rup)); + set_sticky_params( &(currentDirectionalType->w0), + &(currentDirectionalType->v0), + &(currentDirectionalType->v0p), + &(currentDirectionalType->rl), + &(currentDirectionalType->ru), + &(currentDirectionalType->rlp), + &(currentDirectionalType->rup)); } + if( currentAtomType->name[0] != '\0' ){ + isError = 0; + makeAtype( &(currentAtomType->ident), + &(currentAtomType->isLJ), + &(currentDirectionalType->isSticky), + &(currentDirectionalType->isDipole), + &isGB, + &isEAM, + &(currentAtomType->isCharge), + &(currentAtomType->epslon), + &(currentAtomType->sigma), + &(currentAtomType->charge), + &(currentDirectionalType->dipole), + &isError ); + if( isError ){ + sprintf( painCave.errMsg, + "Error initializing the \"%s\" atom type in fortran\n", + currentAtomType->name ); + painCave.isFatal = 1; + simError(); + } + } + currentDirectionalType->next; } - if( currentAtomType->name[0] != '\0' ){ - isError = 0; - makeAtype( &(currentAtomType->ident), - &isGB, - &isEAM, - &(currentAtomType->isDirectional), - &(currentAtomType->isLJ), - &(currentAtomType->isCharge), - &(currentAtomType->epslon), - &(currentAtomType->sigma), - &(currentAtomType->charge), - &isError ); - if( isError ){ - sprintf( painCave.errMsg, - "Error initializing the \"%s\" atom type in fortran\n", - currentAtomType->name ); - painCave.isFatal = 1; - simError(); + + else { + // use all dummy variables if this is not a directional atom + if( currentAtomType->name[0] != '\0' ){ + isError = 0; + makeAtype( &(currentAtomType->ident), + &(currentAtomType->isLJ), + ¬SSD, + ¬Dipole, + &isGB, + &isEAM, + &(currentAtomType->isCharge), + &(currentAtomType->epslon), + &(currentAtomType->sigma), + &(currentAtomType->charge), + &noDipMoment, + &isError ); + if( isError ){ + sprintf( painCave.errMsg, + "Error initializing the \"%s\" atom type in fortran\n", + currentAtomType->name ); + painCave.isFatal = 1; + simError(); + } } } currentAtomType = currentAtomType->next; @@ -587,106 +630,24 @@ void WATER::readParams( void ){ #ifdef IS_MPI sprintf( checkPointMsg, - "WATER atom structures successfully sent to fortran\n" ); + "WATER atom and directional structures successfully" + "sent to fortran\n" ); MPIcheckPoint(); #endif // is_mpi - - - - // now read in the directional stuff - -#ifdef IS_MPI - if( worldRank == 0 ) { -#endif - // read in the directional types - - headDirectionalType = new LinkedDirectionalType; - - fastForward( "DirectionalTypes", "initializeDirectionals" ); - - // we are now at the directionalTypes section - - eof_test = fgets( readLine, sizeof(readLine), frcFile ); - lineNum++; - - - // read a line, and start parsing out the atom types - - if( eof_test == NULL ){ - sprintf( painCave.errMsg, - "Error in reading directionals from force file at line %d.\n", - lineNum ); - painCave.isFatal = 1; - simError(); - } - - // 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( parseDirectional( readLine, lineNum, directionalInfo ) ){ - headDirectionalType->add( directionalInfo ); - } - } - eof_test = fgets( readLine, sizeof(readLine), frcFile ); - lineNum++; - } - -#ifdef IS_MPI - - // send out the linked list to all the other processes - - sprintf( checkPointMsg, - "DUFF directional structures read successfully." ); - MPIcheckPoint(); - - currentDirectionalType = headDirectionalType->next; - while( currentDirectionalType != NULL ){ - currentDirectionalType->duplicate( directionalInfo ); - sendFrcStruct( &directionalInfo, mpiDirectionalStructType ); - currentDirectionalType = currentDirectionalType->next; - } - directionalInfo.last = 1; - sendFrcStruct( &directionalInfo, mpiDirectionalStructType ); - - } - - else{ - - // listen for node 0 to send out the force params - - MPIcheckPoint(); - - headDirectionalType = new LinkedDirectionalType; - receiveFrcStruct( &directionalInfo, mpiDirectionalStructType ); - while( !directionalInfo.last ){ - - headDirectionalType->add( directionalInfo ); - receiveFrcStruct( &directionalInfo, mpiDirectionalStructType ); - } - } - - sprintf( checkPointMsg, - "WATER directional structures broadcast successfully." ); - MPIcheckPoint(); - -#endif // is_mpi } void WATER::initializeAtoms( int nAtoms, Atom** the_atoms ){ - int i; + int i,j,k; // initialize the atoms - + DirectionalAtom* dAtom; + double ji[3]; + double inertialMat[3][3]; for( i=0; ifind( the_atoms[i]->getType() ); if( currentAtomType == NULL ){ sprintf( painCave.errMsg, @@ -695,44 +656,51 @@ void WATER::initializeAtoms( int nAtoms, Atom** the_at painCave.isFatal = 1; simError(); } - - the_atoms[i]->setisLJ( currentAtomType->isLJ); - the_atoms[i]->setisCharge( currentAtomType->isCharge); the_atoms[i]->setMass( currentAtomType->mass ); - the_atoms[i]->setEpslon( currentAtomType->epslon ); - the_atoms[i]->setSigma( currentAtomType->sigma ); - the_atoms[i]->setCharge( currentAtomType->charge ); the_atoms[i]->setIdent( currentAtomType->ident ); - } -} + if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma; -void WATER::initializeDirectional( int nAtoms, Atom** the_atoms ){ - - int i; - - // initialize the atoms - + if( currentAtomType->isDirectional ){ + currentDirectionalType = + headDirectionalType->find( the_atoms[i]->getType() ); + if( currentDirectionalType == NULL ){ + sprintf( painCave.errMsg, + "DirectionalType error, %s not found in force file.\n", + the_atoms[i]->getType() ); + painCave.isFatal = 1; + simError(); + } - for( i=0; ifind( the_atoms[i]->getType() ); - if( currentAtomType == NULL ){ - sprintf( painCave.errMsg, - "AtomType error, %s not found in force file.\n", - the_atoms[i]->getType() ); + // zero out the moments of inertia matrix + for( j=0; j<3; j++ ) + for( k=0; k<3; k++ ) + inertialMat[j][k] = 0.0; + + // load the force file moments of inertia + inertialMat[0][0] = currentDirectionalType->Ixx; + inertialMat[1][1] = currentDirectionalType->Iyy; + inertialMat[2][2] = currentDirectionalType->Izz; + + dAtom = (DirectionalAtom *) the_atoms[i]; + dAtom->setHasDipole( currentDirectionalType->isDipole ); + + ji[0] = 0.0; + ji[1] = 0.0; + ji[2] = 0.0; + dAtom->setJ( ji ); + dAtom->setI( inertialMat ); + + entry_plug->n_dipoles++; + } + else{ + sprintf( painCave.errMsg, + "WATER error: Atom \"%s\" is directional, yet no standard" + " orientation was specifed in the BASS file.\n", + currentAtomType->name ); painCave.isFatal = 1; simError(); } - - the_atoms[i]->setisLJ( currentAtomType->isLJ); - the_atoms[i]->setisCharge( currentAtomType->isCharge); - the_atoms[i]->setMass( currentAtomType->mass ); - the_atoms[i]->setEpslon( currentAtomType->epslon ); - the_atoms[i]->setSigma( currentAtomType->sigma ); - the_atoms[i]->setCharge( currentAtomType->charge ); - the_atoms[i]->setIdent( currentAtomType->ident ); - } } @@ -828,10 +796,11 @@ void WATER::sectionSearch( char* secHead, char* stopTe int foundSection = 0; int foundText = 0; char* the_token; - tempPos = new fpos_t; + fpos_t *tempPos; rewind( frcFile ); lineNum = 0; + tempPos = new fpos_t; eof_test = fgets( readLine, sizeof(readLine), frcFile ); lineNum++; @@ -843,7 +812,6 @@ void WATER::sectionSearch( char* secHead, char* stopTe simError(); } - while( !foundSection ){ while( eof_test != NULL && readLine[0] != '#' ){ eof_test = fgets( readLine, sizeof(readLine), frcFile ); @@ -858,11 +826,9 @@ void WATER::sectionSearch( char* secHead, char* stopTe painCave.isFatal = 1; simError(); } - the_token = strtok( readLine, " ,;\t#\n" ); foundSection = !strcmp( secHead, the_token ); - if( !foundSection ){ eof_test = fgets( readLine, sizeof(readLine), frcFile ); lineNum++; @@ -877,30 +843,32 @@ void WATER::sectionSearch( char* secHead, char* stopTe simError(); } } - } + } while( !foundText ){ - if( !foundText ){ - fgetpos( frcFile, tempPos ); - 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; + + fgetpos( frcFile, tempPos ); + 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( the_token != NULL ){ + foundText = !strcmp( stopText, the_token ); + } } - - fsetPos( frcFile, tempPos ); + fsetpos( frcFile, tempPos ); + eof_test = fgets( readLine, sizeof(readLine), frcFile ); + lineNum++; } @@ -1015,7 +983,7 @@ int WATER_NS::parseDirectional( char *lineBuffer, int simError(); } - info.I_xx = atof( the_token ); + info.Ixx = atof( the_token ); if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ sprintf( painCave.errMsg, @@ -1024,7 +992,7 @@ int WATER_NS::parseDirectional( char *lineBuffer, int simError(); } - info.I_yy = atof( the_token ); + info.Iyy = atof( the_token ); if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ sprintf( painCave.errMsg, @@ -1033,7 +1001,7 @@ int WATER_NS::parseDirectional( char *lineBuffer, int simError(); } - info.I_zz = atof( the_token ); + info.Izz = atof( the_token ); if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ sprintf( painCave.errMsg,