--- trunk/OOPSE/libmdtools/DUFF.cpp 2003/06/19 22:02:44 559 +++ trunk/OOPSE/libmdtools/DUFF.cpp 2003/08/20 19:11:51 704 @@ -15,9 +15,16 @@ using namespace std; #include "mpiForceField.h" #endif // is_mpi -namespace TPE { // restrict the access of the folowing to this file only. +// define some bond Types +#define FIXED_BOND 0 +#define HARMONIC_BOND 1 + + +namespace DUFF_NS { // restrict the access of the folowing to this file only. + + // Declare the structures that will be passed by MPI typedef struct{ @@ -28,6 +35,11 @@ namespace TPE { // restrict the access of the folowin double dipole; double w0; double v0; + double v0p; + double rl; + double ru; + double rlp; + double rup; int isSSD; int isDipole; int ident; @@ -39,10 +51,11 @@ namespace TPE { // restrict the access of the folowin typedef struct{ char nameA[15]; char nameB[15]; - char type[30]; double d0; + double k0; int last; // 0 -> default // 1 -> tells nodes to stop listening + int type; } bondStruct; @@ -111,8 +124,8 @@ namespace TPE { // restrict the access of the folowin if( !strcmp( info.name, name ) ){ sprintf( painCave.errMsg, - "Duplicate TraPPE_Ex atom type \"%s\" found in " - "the TraPPE_ExFF param file./n", + "Duplicate DUFF atom type \"%s\" found in " + "the DUFF param file./n", name ); painCave.isFatal = 1; simError(); @@ -130,6 +143,11 @@ namespace TPE { // restrict the access of the folowin next->dipole = info.dipole; next->w0 = info.w0; next->v0 = info.v0; + next->v0p = info.v0p; + next->rl = info.rl; + next->ru = info.ru; + next->rlp = info.rlp; + next->rup = info.rup; next->ident = info.ident; } } @@ -146,6 +164,11 @@ namespace TPE { // restrict the access of the folowin info.dipole = dipole; info.w0 = w0; info.v0 = v0; + info.v0p = v0p; + info.rl = rl; + info.ru = ru; + info.rlp = rlp; + info.rup = rup; info.ident = ident; info.last = 0; } @@ -162,6 +185,11 @@ namespace TPE { // restrict the access of the folowin double dipole; double w0; double v0; + double v0p; + double rl; + double ru; + double rlp; + double rup; int ident; LinkedAtomType* next; }; @@ -172,7 +200,6 @@ namespace TPE { // restrict the access of the folowin next = NULL; nameA[0] = '\0'; nameB[0] = '\0'; - type[0] = '\0'; } ~LinkedBondType(){ if( next != NULL ) delete next; } @@ -194,8 +221,8 @@ namespace TPE { // restrict the access of the folowin if(dup){ sprintf( painCave.errMsg, - "Duplicate TraPPE_Ex bond type \"%s - %s\" found in " - "the TraPPE_ExFF param file./n", + "Duplicate DUFF bond type \"%s - %s\" found in " + "the DUFF param file./n", nameA, nameB ); painCave.isFatal = 1; simError(); @@ -207,8 +234,9 @@ namespace TPE { // restrict the access of the folowin next = new LinkedBondType(); strcpy(next->nameA, info.nameA); strcpy(next->nameB, info.nameB); - strcpy(next->type, info.type); + next->type = info.type; next->d0 = info.d0; + next->k0 = info.k0; } } @@ -216,8 +244,9 @@ namespace TPE { // restrict the access of the folowin void duplicate( bondStruct &info ){ strcpy(info.nameA, nameA); strcpy(info.nameB, nameB); - strcpy(info.type, type); + info.type = type; info.d0 = d0; + info.k0 = k0; info.last = 0; } @@ -226,8 +255,9 @@ namespace TPE { // restrict the access of the folowin char nameA[15]; char nameB[15]; - char type[30]; + int type; double d0; + double k0; LinkedBondType* next; }; @@ -264,8 +294,8 @@ namespace TPE { // restrict the access of the folowin if(dup){ sprintf( painCave.errMsg, - "Duplicate TraPPE_Ex bend type \"%s - %s - %s\" found in " - "the TraPPE_ExFF param file./n", + "Duplicate DUFF bend type \"%s - %s - %s\" found in " + "the DUFF param file./n", nameA, nameB, nameC ); painCave.isFatal = 1; simError(); @@ -349,8 +379,8 @@ namespace TPE { // restrict the access of the folowin if(dup){ sprintf( painCave.errMsg, - "Duplicate TraPPE_Ex torsion type \"%s - %s - %s - %s\" found in " - "the TraPPE_ExFF param file./n", nameA, nameB, nameC, nameD ); + "Duplicate DUFF torsion type \"%s - %s - %s - %s\" found in " + "the DUFF param file./n", nameA, nameB, nameC, nameD ); painCave.isFatal = 1; simError(); } @@ -410,7 +440,7 @@ namespace TPE { // restrict the access of the folowin } // namespace -using namespace TPE; +using namespace DUFF_NS; //**************************************************************** @@ -418,7 +448,7 @@ using namespace TPE; //**************************************************************** -TraPPE_ExFF::TraPPE_ExFF(){ +DUFF::DUFF(){ char fileName[200]; char* ffPath_env = "FORCE_PARAM_PATH"; @@ -468,7 +498,7 @@ TraPPE_ExFF::TraPPE_ExFF(){ // Init the bondStruct mpi type bondStruct bondProto; // mpiPrototype - int bondBC[3] = {60,1,1}; // block counts + int bondBC[3] = {30,2,2}; // block counts MPI_Aint bondDspls[3]; // displacements MPI_Datatype bondMbrTypes[3]; // member mpi types @@ -537,7 +567,7 @@ TraPPE_ExFF::TraPPE_ExFF(){ // generate the force file name - strcpy( fileName, "TraPPE_Ex.frc" ); + strcpy( fileName, "DUFF.frc" ); // fprintf( stderr,"Trying to open %s\n", fileName ); // attempt to open the file in the current directory first. @@ -576,14 +606,14 @@ TraPPE_ExFF::TraPPE_ExFF(){ #ifdef IS_MPI } - sprintf( checkPointMsg, "TraPPE_ExFF file opened sucessfully." ); + sprintf( checkPointMsg, "DUFF file opened sucessfully." ); MPIcheckPoint(); #endif // is_mpi } -TraPPE_ExFF::~TraPPE_ExFF(){ +DUFF::~DUFF(){ if( headAtomType != NULL ) delete headAtomType; if( headBondType != NULL ) delete headBondType; @@ -601,7 +631,7 @@ TraPPE_ExFF::~TraPPE_ExFF(){ #endif // is_mpi } -void TraPPE_ExFF::cleanMe( void ){ +void DUFF::cleanMe( void ){ #ifdef IS_MPI @@ -620,13 +650,13 @@ void TraPPE_ExFF::cleanMe( void ){ } -void TraPPE_ExFF::initForceField( int ljMixRule ){ +void DUFF::initForceField( int ljMixRule ){ initFortran( ljMixRule, entry_plug->useReactionField ); } -void TraPPE_ExFF::readParams( void ){ +void DUFF::readParams( void ){ int i, a, b, c, d; int identNum; @@ -698,7 +728,7 @@ void TraPPE_ExFF::readParams( void ){ // send out the linked list to all the other processes sprintf( checkPointMsg, - "TraPPE_ExFF atom structures read successfully." ); + "DUFF atom structures read successfully." ); MPIcheckPoint(); currentAtomType = headAtomType->next; //skip the first element who is a place holder. @@ -710,7 +740,7 @@ void TraPPE_ExFF::readParams( void ){ sendFrcStruct( &atomInfo, mpiAtomStructType ); sprintf( checkPointMsg, - "successfully sent TraPPE_Ex force type: \"%s\"\n", + "successfully sent DUFF force type: \"%s\"\n", atomInfo.name ); MPIcheckPoint(); @@ -754,6 +784,7 @@ void TraPPE_ExFF::readParams( void ){ int isGB = 0; int isLJ = 1; + int isEAM =0; double GB_dummy = 0.0; @@ -763,7 +794,10 @@ void TraPPE_ExFF::readParams( void ){ if(currentAtomType->isDipole) entry_plug->useDipole = 1; if(currentAtomType->isSSD) { entry_plug->useSticky = 1; - set_sticky_params( &(currentAtomType->w0), &(currentAtomType->v0)); + set_sticky_params( &(currentAtomType->w0), &(currentAtomType->v0), + &(currentAtomType->v0p), + &(currentAtomType->rl), &(currentAtomType->ru), + &(currentAtomType->rlp), &(currentAtomType->rup)); } if( currentAtomType->name[0] != '\0' ){ @@ -773,6 +807,7 @@ void TraPPE_ExFF::readParams( void ){ &(currentAtomType->isSSD), &(currentAtomType->isDipole), &isGB, + &isEAM, &(currentAtomType->epslon), &(currentAtomType->sigma), &(currentAtomType->dipole), @@ -790,7 +825,7 @@ void TraPPE_ExFF::readParams( void ){ #ifdef IS_MPI sprintf( checkPointMsg, - "TraPPE_ExFF atom structures successfully sent to fortran\n" ); + "DUFF atom structures successfully sent to fortran\n" ); MPIcheckPoint(); #endif // is_mpi @@ -844,7 +879,7 @@ void TraPPE_ExFF::readParams( void ){ // send out the linked list to all the other processes sprintf( checkPointMsg, - "TraPPE_Ex bond structures read successfully." ); + "DUFF bond structures read successfully." ); MPIcheckPoint(); currentBondType = headBondType->next; @@ -874,7 +909,7 @@ void TraPPE_ExFF::readParams( void ){ } sprintf( checkPointMsg, - "TraPPE_ExFF bond structures broadcast successfully." ); + "DUFF bond structures broadcast successfully." ); MPIcheckPoint(); #endif // is_mpi @@ -927,7 +962,7 @@ void TraPPE_ExFF::readParams( void ){ // send out the linked list to all the other processes sprintf( checkPointMsg, - "TraPPE_Ex bend structures read successfully." ); + "DUFF bend structures read successfully." ); MPIcheckPoint(); currentBendType = headBendType->next; @@ -957,7 +992,7 @@ void TraPPE_ExFF::readParams( void ){ } sprintf( checkPointMsg, - "TraPPE_ExFF bend structures broadcast successfully." ); + "DUFF bend structures broadcast successfully." ); MPIcheckPoint(); #endif // is_mpi @@ -1012,7 +1047,7 @@ void TraPPE_ExFF::readParams( void ){ // send out the linked list to all the other processes sprintf( checkPointMsg, - "TraPPE_Ex torsion structures read successfully." ); + "DUFF torsion structures read successfully." ); MPIcheckPoint(); currentTorsionType = headTorsionType->next; @@ -1042,7 +1077,7 @@ void TraPPE_ExFF::readParams( void ){ } sprintf( checkPointMsg, - "TraPPE_ExFF torsion structures broadcast successfully." ); + "DUFF torsion structures broadcast successfully." ); MPIcheckPoint(); #endif // is_mpi @@ -1052,7 +1087,7 @@ void TraPPE_ExFF::readParams( void ){ -void TraPPE_ExFF::initializeAtoms( int nAtoms, Atom** the_atoms ){ +void DUFF::initializeAtoms( int nAtoms, Atom** the_atoms ){ ////////////////////////////////////////////////// @@ -1141,7 +1176,7 @@ void TraPPE_ExFF::initializeAtoms( int nAtoms, Atom** else{ sprintf( painCave.errMsg, - "TraPPE_ExFF error: Atom \"%s\" is a dipole, yet no standard" + "DUFF error: Atom \"%s\" is a dipole, yet no standard" " orientation was specifed in the BASS file.\n", currentAtomType->name ); painCave.isFatal = 1; @@ -1151,7 +1186,7 @@ void TraPPE_ExFF::initializeAtoms( int nAtoms, Atom** else{ if( the_atoms[i]->isDirectional() ){ sprintf( painCave.errMsg, - "TraPPE_ExFF error: Atom \"%s\" was given a standard" + "DUFF error: Atom \"%s\" was given a standard" "orientation in the BASS file, yet it is not a dipole.\n", currentAtomType->name); painCave.isFatal = 1; @@ -1161,7 +1196,7 @@ void TraPPE_ExFF::initializeAtoms( int nAtoms, Atom** } } -void TraPPE_ExFF::initializeBonds( int nBonds, Bond** bondArray, +void DUFF::initializeBonds( int nBonds, Bond** bondArray, bond_pair* the_bonds ){ int i,a,b; char* atomA; @@ -1189,17 +1224,33 @@ void TraPPE_ExFF::initializeBonds( int nBonds, Bond** simError(); } - if( !strcmp( currentBondType->type, "fixed" ) ){ - + switch( currentBondType->type ){ + + case FIXED_BOND: + bondArray[i] = new ConstrainedBond( *the_atoms[a], *the_atoms[b], currentBondType->d0 ); entry_plug->n_constraints++; + break; + + case HARMONIC_BOND: + + bondArray[i] = new HarmonicBond( *the_atoms[a], + *the_atoms[b], + currentBondType->d0, + currentBondType->k0 ); + break; + + default: + + break; + // do nothing } } } -void TraPPE_ExFF::initializeBends( int nBends, Bend** bendArray, +void DUFF::initializeBends( int nBends, Bend** bendArray, bend_set* the_bends ){ QuadraticBend* qBend; @@ -1259,7 +1310,9 @@ void TraPPE_ExFF::initializeBends( int nBends, Bend** } gBend = new GhostBend( *the_atoms[a], - *the_atoms[b] ); + *the_atoms[b], + *the_atoms[c] ); + gBend->setConstants( currentBendType->k1, currentBendType->k2, currentBendType->k3, @@ -1275,12 +1328,12 @@ void TraPPE_ExFF::initializeBends( int nBends, Bend** currentBendType->k3, currentBendType->t0 ); bendArray[i] = qBend; - } + } } } } -void TraPPE_ExFF::initializeTorsions( int nTorsions, Torsion** torsionArray, +void DUFF::initializeTorsions( int nTorsions, Torsion** torsionArray, torsion_set* the_torsions ){ int i, a, b, c, d; @@ -1327,7 +1380,7 @@ void TraPPE_ExFF::initializeTorsions( int nTorsions, T } } -void TraPPE_ExFF::fastForward( char* stopText, char* searchOwner ){ +void DUFF::fastForward( char* stopText, char* searchOwner ){ int foundText = 0; char* the_token; @@ -1382,7 +1435,7 @@ void TraPPE_ExFF::fastForward( char* stopText, char* s } -int TPE::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){ +int DUFF_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){ char* the_token; @@ -1468,17 +1521,62 @@ int TPE::parseAtom( char *lineBuffer, int lineNum, ato } info.v0 = 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.v0p = 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.rl = 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.ru = 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.rlp = 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.rup = atof( the_token ); } - else info.v0 = info.w0 = 0.0; + else info.v0 = info.w0 = info.v0p = info.rl = info.ru = info.rlp = info.rup = 0.0; return 1; } else return 0; } -int TPE::parseBond( char *lineBuffer, int lineNum, bondStruct &info ){ +int DUFF_NS::parseBond( char *lineBuffer, int lineNum, bondStruct &info ){ char* the_token; + char bondType[30]; the_token = strtok( lineBuffer, " \n\t,;" ); if( the_token != NULL ){ @@ -1501,9 +1599,11 @@ int TPE::parseBond( char *lineBuffer, int lineNum, bon simError(); } - strcpy( info.type, the_token ); + strcpy( bondType, the_token ); - if( !strcmp( info.type, "fixed" ) ){ + if( !strcmp( bondType, "fixed" ) ){ + info.type = FIXED_BOND; + if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ sprintf( painCave.errMsg, "Error parseing BondTypes: line %d\n", lineNum ); @@ -1512,10 +1612,34 @@ int TPE::parseBond( char *lineBuffer, int lineNum, bon } info.d0 = atof( the_token ); + + info.k0=0.0; } + else if( !strcmp( bondType, "harmonic" ) ){ + info.type = HARMONIC_BOND; + + if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ + sprintf( painCave.errMsg, + "Error parseing BondTypes: line %d\n", lineNum ); + painCave.isFatal = 1; + simError(); + } + + info.d0 = atof( the_token ); + + if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ + sprintf( painCave.errMsg, + "Error parseing BondTypes: line %d\n", lineNum ); + painCave.isFatal = 1; + simError(); + } + + info.k0 = atof( the_token ); + } + else{ sprintf( painCave.errMsg, - "Unknown TraPPE_Ex bond type \"%s\" at line %d\n", + "Unknown DUFF bond type \"%s\" at line %d\n", info.type, lineNum ); painCave.isFatal = 1; @@ -1528,7 +1652,7 @@ int TPE::parseBond( char *lineBuffer, int lineNum, bon } -int TPE::parseBend( char *lineBuffer, int lineNum, bendStruct &info ){ +int DUFF_NS::parseBend( char *lineBuffer, int lineNum, bendStruct &info ){ char* the_token; @@ -1604,7 +1728,7 @@ int TPE::parseBend( char *lineBuffer, int lineNum, ben else{ sprintf( painCave.errMsg, - "Unknown TraPPE_Ex bend type \"%s\" at line %d\n", + "Unknown DUFF bend type \"%s\" at line %d\n", info.type, lineNum ); painCave.isFatal = 1; @@ -1616,7 +1740,7 @@ int TPE::parseBend( char *lineBuffer, int lineNum, ben else return 0; } -int TPE::parseTorsion( char *lineBuffer, int lineNum, torsionStruct &info ){ +int DUFF_NS::parseTorsion( char *lineBuffer, int lineNum, torsionStruct &info ){ char* the_token; @@ -1702,7 +1826,7 @@ int TPE::parseTorsion( char *lineBuffer, int lineNum, else{ sprintf( painCave.errMsg, - "Unknown TraPPE_Ex torsion type \"%s\" at line %d\n", + "Unknown DUFF torsion type \"%s\" at line %d\n", info.type, lineNum ); painCave.isFatal = 1;