mdtools/interface_implementation/DipoleTestFF.cpp

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

#include <iostream>
using namespace std;

#include "ForceFields.hpp"
#include "SRI.hpp"


DipoleTestFF::DipoleTestFF(){

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