mdtools/libmdCode/DipoleTestFF.cpp

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

#include <iostream>
using namespace std;

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


DipoleTestFF::DipoleTestFF(){

#ifdef IS_MPI
  sprintf( painCave.errMsg,
           "DipoleTest ForceField does not currently support mpi" );
  painCave.isFatal = 1;
  simError();
#endif // is_mpi

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

  // generate the force file name

  strcpy( fileName, "DipoleTest.frc" ); 

  // 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 );
    strcpy( temp, ffPath );
    strcat( temp, "/" );
    strcat( temp, fileName );
    strcpy( fileName, temp );

    frcFile = fopen( fileName, "r" );
    
    if( frcFile == NULL ){
      
      fprintf( stderr,
               "Error opening the force field parameter file: %s\n"
               "Have you tried setting the FORCE_PARAM_PATH environment "
               "vairable?\n",
               fileName );
      exit( 8 );
    }
  }
}

DipoleTestFF::~DipoleTestFF(){
  
  fclose( frcFile );
}

void DipoleTestFF::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;
    }
    
    char name[15];
    int isDipole;
    double mass;
    double epslon;
    double sigma;
    double dipole;
    LinkedType* next;
  };
  
  LinkedType* headAtomType; 
  LinkedType* currentAtomType;
  LinkedType* tempAtomType;

  char readLine[500];
  char* the_token;
  char* eof_test;
  int foundAtom = 0;
  int lineNum = 0;
  int i;
  
  //////////////////////////////////////////////////
  // a quick water fix 

  double testI[3][3]; 
  testI[0][0] = 1.76958347772500;
  testI[0][1] = 0.0;
  testI[0][2] = 0.0;

  testI[1][0] = 0.0;
  testI[1][1] = 0.614537057924513;
  testI[1][2] = 0.0;

  testI[2][0] = 0.0;
  testI[2][1] = 0.0;
  testI[2][2] = 1.15504641980049;

  //////////////////////////////////////////////////

  Atom** the_atoms;
  int nAtoms;
  the_atoms = entry_plug->atoms;
  nAtoms = entry_plug->n_atoms;
  
  // read in the atom types.

  rewind( frcFile );
  headAtomType = new LinkedType;
  currentAtomType = headAtomType;
  
  eof_test = fgets( readLine, sizeof(readLine), frcFile );
  lineNum++;
  if( eof_test == NULL ){
    fprintf( stderr, "Error in reading Atoms from force file.\n" );
    exit(8);
  }

  
  while( !foundAtom ){
    while( eof_test != NULL && readLine[0] != '#' ){
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;
    }
    if( eof_test == NULL ){
      fprintf( stderr, 
               "Error in reading Atoms from force file at line %d.\n",
               lineNum );
      exit(8);
    }

    the_token = strtok( readLine, " ,;\t#\n" );
    foundAtom = !strcmp( "AtomTypes", the_token );
    
    if( !foundAtom ){
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;

      if( eof_test == NULL ){
        fprintf( stderr, 
                 "Error in reading Atoms from force file at line %d.\n",
                 lineNum );
        exit(8);
      } 
    }
  }

  // we are now at the AtomTypes section.

  eof_test = fgets( readLine, sizeof(readLine), frcFile );
  lineNum++;
  
  if( eof_test == NULL ){
    fprintf( stderr, 
             "Error in reading Atoms from force file at line %d.\n",
             lineNum );
    exit(8);
  }

  while( readLine[0] != '#' && eof_test != NULL ){
    
    if( readLine[0] != '!' ){
      
      the_token = strtok( readLine, " \n\t,;" );
      if( the_token != NULL ){
        
        strcpy( currentAtomType->name, the_token );

        if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
          fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
          exit(8);
        }

        sscanf( the_token, "%d", &currentAtomType->isDipole );

        if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
          fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
          exit(8);
        }

        sscanf( the_token, "%lf", &currentAtomType->mass );

        if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
          fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
          exit(8);
        }
        
        sscanf( the_token, "%lf", &currentAtomType->epslon );
        
        if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
          fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
          exit(8);
        }

        sscanf( the_token, "%lf", &currentAtomType->sigma );
        
        if( currentAtomType->isDipole ){
          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
            exit(8);
          }
          
          sscanf( the_token, "%lf", &currentAtomType->dipole );
        }
      }
    }
    
    tempAtomType = new LinkedType;
    currentAtomType->next = tempAtomType;
    currentAtomType = tempAtomType;
    
    eof_test = fgets( readLine, sizeof(readLine), frcFile );
    lineNum++;
  }

  
  // initialize the atoms
  
  DirectionalAtom* dAtom;

  for( i=0; i<nAtoms; i++ ){
    
    currentAtomType = headAtomType->find( the_atoms[i]->getType() );
    if( currentAtomType == NULL ){
      fprintf( stderr, "AtomType error, %s not found in force file.\n",
               the_atoms[i]->getType() );
      exit(8);
    }
    
    the_atoms[i]->setMass( currentAtomType->mass );
    the_atoms[i]->setEpslon( currentAtomType->epslon );
    the_atoms[i]->setSigma( currentAtomType->sigma );
    the_atoms[i]->setLJ();

    if( currentAtomType->isDipole ){
      if( the_atoms[i]->isDirectional() ){
        dAtom = (DirectionalAtom *) the_atoms[i];
        dAtom->setMu( currentAtomType->dipole );
        dAtom->setHasDipole( 1 );
        dAtom->setSSD( 1 );
        dAtom->setJx( 0.0 );
        dAtom->setJy( 0.0 );
        dAtom->setJz( 0.0 );
        
        dAtom->setI( testI );

        entry_plug->n_dipoles++;
      }
      else{
        std::cerr 
          << "DipoleTestFF error: Atom \"" 
          << currentAtomType->name << "\" is a dipole, yet no standard"
          << " orientation was specifed in the BASS file.\n";
        exit(8);
      }
    }
    else{
      if( the_atoms[i]->isDirectional() ){
        std::cerr
          << "DipoleTestFF error: Atom \""
          << currentAtomType->name << "\" was given a standard orientation"
          << " in the BASS file, yet it is not a dipole.\n";
        exit(8);
      }
    }
  }


  // clean up the memory
  
  delete headAtomType;
}

void DipoleTestFF::initializeBonds( bond_pair* the_bonds ){
  
  if( entry_plug->n_bonds ){
    std::cerr << "DipoleTest does not support bonds!\n";
    exit(8);
  }
}

void DipoleTestFF::initializeBends( bend_set* the_bends ){
  
  if( entry_plug->n_bends ){
    std::cerr << "DipoleTest does not support bends!\n";
    exit(8);
  }
}

void DipoleTestFF::initializeTorsions( torsion_set* the_torsions ){

  if( entry_plug->n_torsions ){  
    std::cerr << "DipoleTest does not support torsions!\n";
    exit(8);
  }
}
Revision:	270
Committed:	Fri Feb 14 17:08:46 2003 UTC (22 years, 8 months ago) by mmeineke
File size:	7051 byte(s)
Log Message:	added libmdCode and a couple help scripts
#	User	Rev	Content
1	mmeineke	270	#include <cstdlib>
2			#include <cstdio>
3			#include <cstring>
4
5			#include <iostream>
6			using namespace std;
7
8			#include "ForceFields.hpp"
9			#include "SRI.hpp"
10			#include "simError.h"
11
12
13			DipoleTestFF::DipoleTestFF(){
14
15			#ifdef IS_MPI
16			sprintf( painCave.errMsg,
17			"DipoleTest ForceField does not currently support mpi" );
18			painCave.isFatal = 1;
19			simError();
20			#endif // is_mpi
21
22			char fileName[200];
23			char* ffPath_env = "FORCE_PARAM_PATH";
24			char* ffPath;
25			char temp[200];
26
27			// generate the force file name
28
29			strcpy( fileName, "DipoleTest.frc" );
30
31			// attempt to open the file in the current directory first.
32
33			frcFile = fopen( fileName, "r" );
34
35			if( frcFile == NULL ){
36
37			// next see if the force path enviorment variable is set
38
39			ffPath = getenv( ffPath_env );
40			strcpy( temp, ffPath );
41			strcat( temp, "/" );
42			strcat( temp, fileName );
43			strcpy( fileName, temp );
44
45			frcFile = fopen( fileName, "r" );
46
47			if( frcFile == NULL ){
48
49			fprintf( stderr,
50			"Error opening the force field parameter file: %s\n"
51			"Have you tried setting the FORCE_PARAM_PATH environment "
52			"vairable?\n",
53			fileName );
54			exit( 8 );
55			}
56			}
57			}
58
59			DipoleTestFF::~DipoleTestFF(){
60
61			fclose( frcFile );
62			}
63
64			void DipoleTestFF::initializeAtoms( void ){
65
66			class LinkedType {
67			public:
68			LinkedType(){
69			next = NULL;
70			name[0] = '\0';
71			}
72			~LinkedType(){ if( next != NULL ) delete next; }
73
74			LinkedType* find(char* key){
75			if( !strcmp(name, key) ) return this;
76			if( next != NULL ) return next->find(key);
77			return NULL;
78			}
79
80			char name[15];
81			int isDipole;
82			double mass;
83			double epslon;
84			double sigma;
85			double dipole;
86			LinkedType* next;
87			};
88
89			LinkedType* headAtomType;
90			LinkedType* currentAtomType;
91			LinkedType* tempAtomType;
92
93			char readLine[500];
94			char* the_token;
95			char* eof_test;
96			int foundAtom = 0;
97			int lineNum = 0;
98			int i;
99
100			//////////////////////////////////////////////////
101			// a quick water fix
102
103			double testI[3][3];
104			testI[0][0] = 1.76958347772500;
105			testI[0][1] = 0.0;
106			testI[0][2] = 0.0;
107
108			testI[1][0] = 0.0;
109			testI[1][1] = 0.614537057924513;
110			testI[1][2] = 0.0;
111
112			testI[2][0] = 0.0;
113			testI[2][1] = 0.0;
114			testI[2][2] = 1.15504641980049;
115
116			//////////////////////////////////////////////////
117
118			Atom** the_atoms;
119			int nAtoms;
120			the_atoms = entry_plug->atoms;
121			nAtoms = entry_plug->n_atoms;
122
123			// read in the atom types.
124
125			rewind( frcFile );
126			headAtomType = new LinkedType;
127			currentAtomType = headAtomType;
128
129			eof_test = fgets( readLine, sizeof(readLine), frcFile );
130			lineNum++;
131			if( eof_test == NULL ){
132			fprintf( stderr, "Error in reading Atoms from force file.\n" );
133			exit(8);
134			}
135
136
137			while( !foundAtom ){
138			while( eof_test != NULL && readLine[0] != '#' ){
139			eof_test = fgets( readLine, sizeof(readLine), frcFile );
140			lineNum++;
141			}
142			if( eof_test == NULL ){
143			fprintf( stderr,
144			"Error in reading Atoms from force file at line %d.\n",
145			lineNum );
146			exit(8);
147			}
148
149			the_token = strtok( readLine, " ,;\t#\n" );
150			foundAtom = !strcmp( "AtomTypes", the_token );
151
152			if( !foundAtom ){
153			eof_test = fgets( readLine, sizeof(readLine), frcFile );
154			lineNum++;
155
156			if( eof_test == NULL ){
157			fprintf( stderr,
158			"Error in reading Atoms from force file at line %d.\n",
159			lineNum );
160			exit(8);
161			}
162			}
163			}
164
165			// we are now at the AtomTypes section.
166
167			eof_test = fgets( readLine, sizeof(readLine), frcFile );
168			lineNum++;
169
170			if( eof_test == NULL ){
171			fprintf( stderr,
172			"Error in reading Atoms from force file at line %d.\n",
173			lineNum );
174			exit(8);
175			}
176
177			while( readLine[0] != '#' && eof_test != NULL ){
178
179			if( readLine[0] != '!' ){
180
181			the_token = strtok( readLine, " \n\t,;" );
182			if( the_token != NULL ){
183
184			strcpy( currentAtomType->name, the_token );
185
186			if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
187			fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
188			exit(8);
189			}
190
191			sscanf( the_token, "%d", &currentAtomType->isDipole );
192
193			if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
194			fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
195			exit(8);
196			}
197
198			sscanf( the_token, "%lf", &currentAtomType->mass );
199
200			if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
201			fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
202			exit(8);
203			}
204
205			sscanf( the_token, "%lf", &currentAtomType->epslon );
206
207			if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
208			fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
209			exit(8);
210			}
211
212			sscanf( the_token, "%lf", &currentAtomType->sigma );
213
214			if( currentAtomType->isDipole ){
215			if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
216			fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
217			exit(8);
218			}
219
220			sscanf( the_token, "%lf", &currentAtomType->dipole );
221			}
222			}
223			}
224
225			tempAtomType = new LinkedType;
226			currentAtomType->next = tempAtomType;
227			currentAtomType = tempAtomType;
228
229			eof_test = fgets( readLine, sizeof(readLine), frcFile );
230			lineNum++;
231			}
232
233
234			// initialize the atoms
235
236			DirectionalAtom* dAtom;
237
238			for( i=0; i<nAtoms; i++ ){
239
240			currentAtomType = headAtomType->find( the_atoms[i]->getType() );
241			if( currentAtomType == NULL ){
242			fprintf( stderr, "AtomType error, %s not found in force file.\n",
243			the_atoms[i]->getType() );
244			exit(8);
245			}
246
247			the_atoms[i]->setMass( currentAtomType->mass );
248			the_atoms[i]->setEpslon( currentAtomType->epslon );
249			the_atoms[i]->setSigma( currentAtomType->sigma );
250			the_atoms[i]->setLJ();
251
252			if( currentAtomType->isDipole ){
253			if( the_atoms[i]->isDirectional() ){
254			dAtom = (DirectionalAtom *) the_atoms[i];
255			dAtom->setMu( currentAtomType->dipole );
256			dAtom->setHasDipole( 1 );
257			dAtom->setSSD( 1 );
258			dAtom->setJx( 0.0 );
259			dAtom->setJy( 0.0 );
260			dAtom->setJz( 0.0 );
261
262			dAtom->setI( testI );
263
264			entry_plug->n_dipoles++;
265			}
266			else{
267			std::cerr
268			<< "DipoleTestFF error: Atom \""
269			<< currentAtomType->name << "\" is a dipole, yet no standard"
270			<< " orientation was specifed in the BASS file.\n";
271			exit(8);
272			}
273			}
274			else{
275			if( the_atoms[i]->isDirectional() ){
276			std::cerr
277			<< "DipoleTestFF error: Atom \""
278			<< currentAtomType->name << "\" was given a standard orientation"
279			<< " in the BASS file, yet it is not a dipole.\n";
280			exit(8);
281			}
282			}
283			}
284
285
286			// clean up the memory
287
288			delete headAtomType;
289			}
290
291			void DipoleTestFF::initializeBonds( bond_pair* the_bonds ){
292
293			if( entry_plug->n_bonds ){
294			std::cerr << "DipoleTest does not support bonds!\n";
295			exit(8);
296			}
297			}
298
299			void DipoleTestFF::initializeBends( bend_set* the_bends ){
300
301			if( entry_plug->n_bends ){
302			std::cerr << "DipoleTest does not support bends!\n";
303			exit(8);
304			}
305			}
306
307			void DipoleTestFF::initializeTorsions( torsion_set* the_torsions ){
308
309			if( entry_plug->n_torsions ){
310			std::cerr << "DipoleTest does not support torsions!\n";
311			exit(8);
312			}
313			}