--- trunk/OOPSE/libmdtools/DUFF.cpp 2003/06/19 22:02:44 559 +++ trunk/OOPSE/libmdtools/DUFF.cpp 2004/01/16 15:01:14 950 @@ -1,6 +1,6 @@ -#include -#include -#include +#include +#include +#include #include using namespace std; @@ -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{ @@ -25,10 +32,17 @@ namespace TPE { // restrict the access of the folowin double mass; double epslon; double sigma; + double charge; double dipole; double w0; double v0; + double v0p; + double rl; + double ru; + double rlp; + double rup; int isSSD; + int isCharge; int isDipole; int ident; int last; // 0 -> default @@ -39,10 +53,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 +126,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 +145,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 +166,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 +187,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 +202,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 +223,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 +236,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 +246,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 +257,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 +296,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 +381,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 +442,7 @@ namespace TPE { // restrict the access of the folowin } // namespace -using namespace TPE; +using namespace DUFF_NS; //**************************************************************** @@ -418,13 +450,12 @@ using namespace TPE; //**************************************************************** -TraPPE_ExFF::TraPPE_ExFF(){ +DUFF::DUFF(){ char fileName[200]; char* ffPath_env = "FORCE_PARAM_PATH"; char* ffPath; char temp[200]; - char errMsg[1000]; headAtomType = NULL; currentAtomType = NULL; @@ -446,7 +477,7 @@ TraPPE_ExFF::TraPPE_ExFF(){ // Init the atomStruct mpi type atomStruct atomProto; // mpiPrototype - int atomBC[3] = {15,6,4}; // block counts + int atomBC[3] = {15,12,5}; // block counts MPI_Aint atomDspls[3]; // displacements MPI_Datatype atomMbrTypes[3]; // member mpi types @@ -468,7 +499,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 +568,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 +607,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 +632,7 @@ TraPPE_ExFF::~TraPPE_ExFF(){ #endif // is_mpi } -void TraPPE_ExFF::cleanMe( void ){ +void DUFF::cleanMe( void ){ #ifdef IS_MPI @@ -620,20 +651,15 @@ 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; - char* atomA; - char* atomB; - char* atomC; - char* atomD; atomStruct atomInfo; bondStruct bondInfo; @@ -698,19 +724,17 @@ 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. while( currentAtomType != NULL ){ currentAtomType->duplicate( atomInfo ); - - sendFrcStruct( &atomInfo, mpiAtomStructType ); sprintf( checkPointMsg, - "successfully sent TraPPE_Ex force type: \"%s\"\n", + "successfully sent DUFF force type: \"%s\"\n", atomInfo.name ); MPIcheckPoint(); @@ -724,7 +748,7 @@ void TraPPE_ExFF::readParams( void ){ else{ // listen for node 0 to send out the force params - + MPIcheckPoint(); headAtomType = new LinkedAtomType; @@ -732,8 +756,6 @@ void TraPPE_ExFF::readParams( void ){ while( !atomInfo.last ){ - - headAtomType->add( atomInfo ); MPIcheckPoint(); @@ -754,16 +776,20 @@ void TraPPE_ExFF::readParams( void ){ int isGB = 0; int isLJ = 1; - double GB_dummy = 0.0; - - + int isEAM =0; + int isCharge = 0; + double charge=0.0; + currentAtomType = headAtomType->next;; while( currentAtomType != NULL ){ - if(currentAtomType->isDipole) entry_plug->useDipole = 1; + if(currentAtomType->isDipole) entry_plug->useDipoles = 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,8 +799,11 @@ void TraPPE_ExFF::readParams( void ){ &(currentAtomType->isSSD), &(currentAtomType->isDipole), &isGB, + &isEAM, + &isCharge, &(currentAtomType->epslon), &(currentAtomType->sigma), + &charge, &(currentAtomType->dipole), &isError ); if( isError ){ @@ -790,7 +819,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 +873,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 +903,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 +956,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 +986,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 +1041,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 +1071,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 +1081,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 +1170,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 +1180,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 +1190,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 +1218,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 +1304,8 @@ void TraPPE_ExFF::initializeBends( int nBends, Bend** } gBend = new GhostBend( *the_atoms[a], - *the_atoms[b] ); + *the_atoms[b]); + gBend->setConstants( currentBendType->k1, currentBendType->k2, currentBendType->k3, @@ -1275,12 +1321,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 +1373,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 +1428,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 +1514,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 +1592,25 @@ 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 ); + painCave.isFatal = 1; + simError(); + } + + 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 ); @@ -1512,11 +1619,21 @@ int TPE::parseBond( char *lineBuffer, int lineNum, bon } 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", - info.type, + "Unknown DUFF bond type \"%s\" at line %d\n", + bondType, lineNum ); painCave.isFatal = 1; simError(); @@ -1528,7 +1645,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 +1721,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 +1733,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 +1819,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;