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 |
TraPPEFF::TraPPEFF(){ |
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, "TraPPE.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 |
TraPPEFF::~TraPPEFF(){ |
52 |
|
53 |
fclose( frcFile ); |
54 |
} |
55 |
|
56 |
void TraPPEFF::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 |
double mass; |
74 |
double epslon; |
75 |
double sigma; |
76 |
LinkedType* next; |
77 |
}; |
78 |
|
79 |
LinkedType* headAtomType; |
80 |
LinkedType* currentAtomType; |
81 |
LinkedType* tempAtomType; |
82 |
|
83 |
char readLine[500]; |
84 |
char* the_token; |
85 |
char* eof_test; |
86 |
int foundAtom = 0; |
87 |
int lineNum = 0; |
88 |
int i; |
89 |
|
90 |
Atom** the_atoms; |
91 |
int nAtoms; |
92 |
the_atoms = entry_plug->atoms; |
93 |
nAtoms = entry_plug->n_atoms; |
94 |
|
95 |
// read in the atom types. |
96 |
|
97 |
rewind( frcFile ); |
98 |
headAtomType = new LinkedType; |
99 |
currentAtomType = headAtomType; |
100 |
|
101 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
102 |
lineNum++; |
103 |
if( eof_test == NULL ){ |
104 |
fprintf( stderr, "Error in reading Atoms from force file.\n" ); |
105 |
exit(8); |
106 |
} |
107 |
|
108 |
|
109 |
while( !foundAtom ){ |
110 |
while( eof_test != NULL && readLine[0] != '#' ){ |
111 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
112 |
lineNum++; |
113 |
} |
114 |
if( eof_test == NULL ){ |
115 |
fprintf( stderr, |
116 |
"Error in reading Atoms from force file at line %d.\n", |
117 |
lineNum ); |
118 |
exit(8); |
119 |
} |
120 |
|
121 |
the_token = strtok( readLine, " ,;\t#\n" ); |
122 |
foundAtom = !strcmp( "AtomTypes", the_token ); |
123 |
|
124 |
if( !foundAtom ){ |
125 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
126 |
lineNum++; |
127 |
|
128 |
if( eof_test == NULL ){ |
129 |
fprintf( stderr, |
130 |
"Error in reading Atoms from force file at line %d.\n", |
131 |
lineNum ); |
132 |
exit(8); |
133 |
} |
134 |
} |
135 |
} |
136 |
|
137 |
// we are now at the AtomTypes section. |
138 |
|
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 |
while( readLine[0] != '#' && eof_test != NULL ){ |
150 |
|
151 |
if( readLine[0] != '!' ){ |
152 |
|
153 |
the_token = strtok( readLine, " \n\t,;" ); |
154 |
if( the_token != NULL ){ |
155 |
|
156 |
strcpy( currentAtomType->name, the_token ); |
157 |
|
158 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
159 |
fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum ); |
160 |
exit(8); |
161 |
} |
162 |
|
163 |
sscanf( the_token, "%lf", ¤tAtomType->mass ); |
164 |
|
165 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
166 |
fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum ); |
167 |
exit(8); |
168 |
} |
169 |
|
170 |
sscanf( the_token, "%lf", ¤tAtomType->epslon ); |
171 |
|
172 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
173 |
fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum ); |
174 |
exit(8); |
175 |
} |
176 |
|
177 |
sscanf( the_token, "%lf", ¤tAtomType->sigma ); |
178 |
} |
179 |
} |
180 |
|
181 |
tempAtomType = new LinkedType; |
182 |
currentAtomType->next = tempAtomType; |
183 |
currentAtomType = tempAtomType; |
184 |
|
185 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
186 |
lineNum++; |
187 |
} |
188 |
|
189 |
|
190 |
// initialize the atoms |
191 |
|
192 |
for( i=0; i<nAtoms; i++ ){ |
193 |
|
194 |
currentAtomType = headAtomType->find( the_atoms[i]->getType() ); |
195 |
if( currentAtomType == NULL ){ |
196 |
fprintf( stderr, "AtomType error, %s not found in force file.\n", |
197 |
the_atoms[i]->getType() ); |
198 |
exit(8); |
199 |
} |
200 |
|
201 |
the_atoms[i]->setMass( currentAtomType->mass ); |
202 |
the_atoms[i]->setEpslon( currentAtomType->epslon ); |
203 |
the_atoms[i]->setSigma( currentAtomType->sigma ); |
204 |
the_atoms[i]->setLJ(); |
205 |
} |
206 |
|
207 |
|
208 |
// clean up the memory |
209 |
|
210 |
delete headAtomType; |
211 |
} |
212 |
|
213 |
void TraPPEFF::initializeBonds( bond_pair* the_bonds ){ |
214 |
|
215 |
class LinkedType { |
216 |
public: |
217 |
LinkedType(){ |
218 |
next = NULL; |
219 |
nameA[0] = '\0'; |
220 |
nameB[0] = '\0'; |
221 |
type[0] = '\0'; |
222 |
} |
223 |
~LinkedType(){ if( next != NULL ) delete next; } |
224 |
|
225 |
LinkedType* find(char* key1, char* key2){ |
226 |
if( !strcmp(nameA, key1 ) && !strcmp( nameB, key2 ) ) return this; |
227 |
if( !strcmp(nameA, key2 ) && !strcmp( nameB, key1 ) ) return this; |
228 |
if( next != NULL ) return next->find(key1, key2); |
229 |
return NULL; |
230 |
} |
231 |
|
232 |
char nameA[15]; |
233 |
char nameB[15]; |
234 |
char type[30]; |
235 |
double d0; |
236 |
|
237 |
LinkedType* next; |
238 |
}; |
239 |
|
240 |
|
241 |
|
242 |
LinkedType* headBondType; |
243 |
LinkedType* currentBondType; |
244 |
LinkedType* tempBondType; |
245 |
|
246 |
char readLine[500]; |
247 |
char* the_token; |
248 |
char* eof_test; |
249 |
int foundBond = 0; |
250 |
int lineNum = 0; |
251 |
int i, a, b; |
252 |
char* atomA; |
253 |
char* atomB; |
254 |
|
255 |
SRI **the_sris; |
256 |
Atom** the_atoms; |
257 |
int nBonds; |
258 |
the_sris = entry_plug->sr_interactions; |
259 |
the_atoms = entry_plug->atoms; |
260 |
nBonds = entry_plug->n_bonds; |
261 |
|
262 |
// read in the bond types. |
263 |
|
264 |
rewind( frcFile ); |
265 |
headBondType = new LinkedType; |
266 |
currentBondType = headBondType; |
267 |
|
268 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
269 |
lineNum++; |
270 |
if( eof_test == NULL ){ |
271 |
fprintf( stderr, "Error in reading Bonds from force file.\n" ); |
272 |
exit(8); |
273 |
} |
274 |
|
275 |
|
276 |
while( !foundBond ){ |
277 |
while( eof_test != NULL && readLine[0] != '#' ){ |
278 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
279 |
lineNum++; |
280 |
} |
281 |
if( eof_test == NULL ){ |
282 |
fprintf( stderr, |
283 |
"Error in reading Bonds from force file at line %d.\n", |
284 |
lineNum ); |
285 |
exit(8); |
286 |
} |
287 |
|
288 |
the_token = strtok( readLine, " ,;\t#\n" ); |
289 |
foundBond = !strcmp( "BondTypes", the_token ); |
290 |
|
291 |
if( !foundBond ){ |
292 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
293 |
lineNum++; |
294 |
|
295 |
if( eof_test == NULL ){ |
296 |
fprintf( stderr, |
297 |
"Error in reading Bonds from force file at line %d.\n", |
298 |
lineNum ); |
299 |
exit(8); |
300 |
} |
301 |
} |
302 |
} |
303 |
|
304 |
// we are now at the BondTypes section. |
305 |
|
306 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
307 |
lineNum++; |
308 |
|
309 |
if( eof_test == NULL ){ |
310 |
fprintf( stderr, |
311 |
"Error in reading Bonds from force file at line %d.\n", |
312 |
lineNum ); |
313 |
exit(8); |
314 |
} |
315 |
|
316 |
while( readLine[0] != '#' && eof_test != NULL ){ |
317 |
|
318 |
if( readLine[0] != '!' ){ |
319 |
|
320 |
the_token = strtok( readLine, " \n\t,;" ); |
321 |
if( the_token != NULL ){ |
322 |
|
323 |
strcpy( currentBondType->nameA, the_token ); |
324 |
|
325 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
326 |
fprintf( stderr, "Error parseing BondTypes: line %d\n", lineNum ); |
327 |
exit(8); |
328 |
} |
329 |
|
330 |
strcpy( currentBondType->nameB, the_token ); |
331 |
|
332 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
333 |
fprintf( stderr, "Error parseing BondTypes: line %d\n", lineNum ); |
334 |
exit(8); |
335 |
} |
336 |
|
337 |
strcpy( currentBondType->type, the_token ); |
338 |
|
339 |
if( !strcmp( currentBondType->type, "fixed" ) ){ |
340 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
341 |
fprintf( stderr, "Error parseing BondTypes: line %d\n", lineNum ); |
342 |
exit(8); |
343 |
} |
344 |
|
345 |
sscanf( the_token, "%lf", ¤tBondType->d0 ); |
346 |
} |
347 |
|
348 |
else{ |
349 |
fprintf(stderr, |
350 |
"Unknown TraPPE bond type \"%s\" at line %d\n", |
351 |
currentBondType->type, |
352 |
lineNum ); |
353 |
exit(8); |
354 |
} |
355 |
} |
356 |
} |
357 |
|
358 |
tempBondType = new LinkedType; |
359 |
currentBondType->next = tempBondType; |
360 |
currentBondType = tempBondType; |
361 |
|
362 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
363 |
lineNum++; |
364 |
} |
365 |
|
366 |
|
367 |
// initialize the Bonds |
368 |
|
369 |
|
370 |
for( i=0; i<nBonds; i++ ){ |
371 |
|
372 |
a = the_bonds[i].a; |
373 |
b = the_bonds[i].b; |
374 |
|
375 |
atomA = the_atoms[a]->getType(); |
376 |
atomB = the_atoms[b]->getType(); |
377 |
currentBondType = headBondType->find( atomA, atomB ); |
378 |
if( currentBondType == NULL ){ |
379 |
fprintf( stderr, "BondType error, %s - %s not found in force file.\n", |
380 |
atomA, atomB ); |
381 |
exit(8); |
382 |
} |
383 |
|
384 |
if( !strcmp( currentBondType->type, "fixed" ) ){ |
385 |
|
386 |
the_sris[i] = new ConstrainedBond( *the_atoms[a], |
387 |
*the_atoms[b], |
388 |
currentBondType->d0 ); |
389 |
entry_plug->n_constraints++; |
390 |
} |
391 |
} |
392 |
|
393 |
|
394 |
// clean up the memory |
395 |
|
396 |
delete headBondType; |
397 |
} |
398 |
|
399 |
void TraPPEFF::initializeBends( bend_set* the_bends ){ |
400 |
|
401 |
class LinkedType { |
402 |
public: |
403 |
LinkedType(){ |
404 |
next = NULL; |
405 |
nameA[0] = '\0'; |
406 |
nameB[0] = '\0'; |
407 |
nameC[0] = '\0'; |
408 |
type[0] = '\0'; |
409 |
} |
410 |
~LinkedType(){ if( next != NULL ) delete next; } |
411 |
|
412 |
LinkedType* find( char* key1, char* key2, char* key3 ){ |
413 |
if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) |
414 |
&& !strcmp( nameC, key3 ) ) return this; |
415 |
if( !strcmp( nameA, key3 ) && !strcmp( nameB, key2 ) |
416 |
&& !strcmp( nameC, key1 ) ) return this; |
417 |
if( next != NULL ) return next->find(key1, key2, key3); |
418 |
return NULL; |
419 |
} |
420 |
|
421 |
char nameA[15]; |
422 |
char nameB[15]; |
423 |
char nameC[15]; |
424 |
char type[30]; |
425 |
double k1, k2, k3, t0; |
426 |
|
427 |
LinkedType* next; |
428 |
}; |
429 |
|
430 |
LinkedType* headBendType; |
431 |
LinkedType* currentBendType; |
432 |
LinkedType* tempBendType; |
433 |
|
434 |
char readLine[500]; |
435 |
char* the_token; |
436 |
char* eof_test; |
437 |
int foundBend = 0; |
438 |
int lineNum = 0; |
439 |
int i, a, b, c, index; |
440 |
char* atomA; |
441 |
char* atomB; |
442 |
char* atomC; |
443 |
QuadraticBend* qBend; |
444 |
|
445 |
SRI **the_sris; |
446 |
Atom** the_atoms; |
447 |
int nBends; |
448 |
the_sris = entry_plug->sr_interactions; |
449 |
the_atoms = entry_plug->atoms; |
450 |
nBends = entry_plug->n_bends; |
451 |
|
452 |
// read in the bend types. |
453 |
|
454 |
rewind( frcFile ); |
455 |
headBendType = new LinkedType; |
456 |
currentBendType = headBendType; |
457 |
|
458 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
459 |
lineNum++; |
460 |
if( eof_test == NULL ){ |
461 |
fprintf( stderr, "Error in reading Bends from force file.\n" ); |
462 |
exit(8); |
463 |
} |
464 |
|
465 |
|
466 |
while( !foundBend ){ |
467 |
while( eof_test != NULL && readLine[0] != '#' ){ |
468 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
469 |
lineNum++; |
470 |
} |
471 |
if( eof_test == NULL ){ |
472 |
fprintf( stderr, |
473 |
"Error in reading Bends from force file at line %d.\n", |
474 |
lineNum ); |
475 |
exit(8); |
476 |
} |
477 |
|
478 |
the_token = strtok( readLine, " ,;\t#\n" ); |
479 |
foundBend = !strcmp( "BendTypes", the_token ); |
480 |
|
481 |
if( !foundBend ){ |
482 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
483 |
lineNum++; |
484 |
|
485 |
if( eof_test == NULL ){ |
486 |
fprintf( stderr, |
487 |
"Error in reading Bends from force file at line %d.\n", |
488 |
lineNum ); |
489 |
exit(8); |
490 |
} |
491 |
} |
492 |
} |
493 |
|
494 |
// we are now at the BendTypes section. |
495 |
|
496 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
497 |
lineNum++; |
498 |
|
499 |
if( eof_test == NULL ){ |
500 |
fprintf( stderr, |
501 |
"Error in reading Bends from force file at line %d.\n", |
502 |
lineNum ); |
503 |
exit(8); |
504 |
} |
505 |
|
506 |
while( readLine[0] != '#' && eof_test != NULL ){ |
507 |
|
508 |
if( readLine[0] != '!' ){ |
509 |
|
510 |
the_token = strtok( readLine, " \n\t,;" ); |
511 |
if( the_token != NULL ){ |
512 |
|
513 |
strcpy( currentBendType->nameA, the_token ); |
514 |
|
515 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
516 |
fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum ); |
517 |
exit(8); |
518 |
} |
519 |
|
520 |
strcpy( currentBendType->nameB, the_token ); |
521 |
|
522 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
523 |
fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum ); |
524 |
exit(8); |
525 |
} |
526 |
|
527 |
strcpy( currentBendType->nameC, the_token ); |
528 |
|
529 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
530 |
fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum ); |
531 |
exit(8); |
532 |
} |
533 |
|
534 |
strcpy( currentBendType->type, the_token ); |
535 |
|
536 |
if( !strcmp( currentBendType->type, "quadratic" ) ){ |
537 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
538 |
fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum ); |
539 |
exit(8); |
540 |
} |
541 |
|
542 |
sscanf( the_token, "%lf", ¤tBendType->k1 ); |
543 |
|
544 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
545 |
fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum ); |
546 |
exit(8); |
547 |
} |
548 |
|
549 |
sscanf( the_token, "%lf", ¤tBendType->k2 ); |
550 |
|
551 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
552 |
fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum ); |
553 |
exit(8); |
554 |
} |
555 |
|
556 |
sscanf( the_token, "%lf", ¤tBendType->k3 ); |
557 |
|
558 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
559 |
fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum ); |
560 |
exit(8); |
561 |
} |
562 |
|
563 |
sscanf( the_token, "%lf", ¤tBendType->t0 ); |
564 |
} |
565 |
|
566 |
else{ |
567 |
fprintf(stderr, |
568 |
"Unknown TraPPE bend type \"%s\" at line %d\n", |
569 |
currentBendType->type, |
570 |
lineNum ); |
571 |
exit(8); |
572 |
} |
573 |
} |
574 |
} |
575 |
|
576 |
tempBendType = new LinkedType; |
577 |
currentBendType->next = tempBendType; |
578 |
currentBendType = tempBendType; |
579 |
|
580 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
581 |
lineNum++; |
582 |
} |
583 |
|
584 |
|
585 |
// initialize the Bends |
586 |
|
587 |
for( i=0; i<nBends; i++ ){ |
588 |
|
589 |
a = the_bends[i].a; |
590 |
b = the_bends[i].b; |
591 |
c = the_bends[i].c; |
592 |
|
593 |
atomA = the_atoms[a]->getType(); |
594 |
atomB = the_atoms[b]->getType(); |
595 |
atomC = the_atoms[c]->getType(); |
596 |
currentBendType = headBendType->find( atomA, atomB, atomC ); |
597 |
if( currentBendType == NULL ){ |
598 |
fprintf( stderr, "BendType error, %s - %s - %s not found" |
599 |
" in force file.\n", |
600 |
atomA, atomB, atomC ); |
601 |
exit(8); |
602 |
} |
603 |
|
604 |
if( !strcmp( currentBendType->type, "quadratic" ) ){ |
605 |
|
606 |
index = i + entry_plug->n_bonds; |
607 |
qBend = new QuadraticBend( *the_atoms[a], |
608 |
*the_atoms[b], |
609 |
*the_atoms[c] ); |
610 |
qBend->setConstants( currentBendType->k1, |
611 |
currentBendType->k2, |
612 |
currentBendType->k3, |
613 |
currentBendType->t0 ); |
614 |
the_sris[index] = qBend; |
615 |
} |
616 |
} |
617 |
|
618 |
|
619 |
// clean up the memory |
620 |
|
621 |
delete headBendType; |
622 |
} |
623 |
|
624 |
void TraPPEFF::initializeTorsions( torsion_set* the_torsions ){ |
625 |
|
626 |
class LinkedType { |
627 |
public: |
628 |
LinkedType(){ |
629 |
next = NULL; |
630 |
nameA[0] = '\0'; |
631 |
nameB[0] = '\0'; |
632 |
nameC[0] = '\0'; |
633 |
type[0] = '\0'; |
634 |
} |
635 |
~LinkedType(){ if( next != NULL ) delete next; } |
636 |
|
637 |
LinkedType* find( char* key1, char* key2, char* key3, char* key4 ){ |
638 |
if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) && |
639 |
!strcmp( nameC, key3 ) && !strcmp( nameD, key4 ) ) return this; |
640 |
|
641 |
if( !strcmp( nameA, key4 ) && !strcmp( nameB, key3 ) && |
642 |
!strcmp( nameC, key2 ) && !strcmp( nameD, key4 ) ) return this; |
643 |
|
644 |
if( next != NULL ) return next->find(key1, key2, key3, key4); |
645 |
return NULL; |
646 |
} |
647 |
|
648 |
char nameA[15]; |
649 |
char nameB[15]; |
650 |
char nameC[15]; |
651 |
char nameD[15]; |
652 |
char type[30]; |
653 |
double k1, k2, k3, k4; |
654 |
|
655 |
LinkedType* next; |
656 |
}; |
657 |
|
658 |
LinkedType* headTorsionType; |
659 |
LinkedType* currentTorsionType; |
660 |
LinkedType* tempTorsionType; |
661 |
|
662 |
char readLine[500]; |
663 |
char* the_token; |
664 |
char* eof_test; |
665 |
int foundTorsion = 0; |
666 |
int lineNum = 0; |
667 |
int i, a, b, c, d, index; |
668 |
char* atomA; |
669 |
char* atomB; |
670 |
char* atomC; |
671 |
char* atomD; |
672 |
CubicTorsion* cTors; |
673 |
|
674 |
SRI **the_sris; |
675 |
Atom** the_atoms; |
676 |
int nTorsions; |
677 |
the_sris = entry_plug->sr_interactions; |
678 |
the_atoms = entry_plug->atoms; |
679 |
nTorsions = entry_plug->n_torsions; |
680 |
|
681 |
// read in the torsion types. |
682 |
|
683 |
rewind( frcFile ); |
684 |
headTorsionType = new LinkedType; |
685 |
currentTorsionType = headTorsionType; |
686 |
|
687 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
688 |
lineNum++; |
689 |
if( eof_test == NULL ){ |
690 |
fprintf( stderr, "Error in reading Torsions from force file.\n" ); |
691 |
exit(8); |
692 |
} |
693 |
|
694 |
|
695 |
while( !foundTorsion ){ |
696 |
while( eof_test != NULL && readLine[0] != '#' ){ |
697 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
698 |
lineNum++; |
699 |
} |
700 |
if( eof_test == NULL ){ |
701 |
fprintf( stderr, |
702 |
"Error in reading Torsions from force file at line %d.\n", |
703 |
lineNum ); |
704 |
exit(8); |
705 |
} |
706 |
|
707 |
the_token = strtok( readLine, " ,;\t#\n" ); |
708 |
foundTorsion = !strcmp( "TorsionTypes", the_token ); |
709 |
|
710 |
if( !foundTorsion ){ |
711 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
712 |
lineNum++; |
713 |
|
714 |
if( eof_test == NULL ){ |
715 |
fprintf( stderr, |
716 |
"Error in reading Torsions from force file at line %d.\n", |
717 |
lineNum ); |
718 |
exit(8); |
719 |
} |
720 |
} |
721 |
} |
722 |
|
723 |
// we are now at the TorsionTypes section. |
724 |
|
725 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
726 |
lineNum++; |
727 |
|
728 |
if( eof_test == NULL ){ |
729 |
fprintf( stderr, |
730 |
"Error in reading Torsions from force file at line %d.\n", |
731 |
lineNum ); |
732 |
exit(8); |
733 |
} |
734 |
|
735 |
while( readLine[0] != '#' && eof_test != NULL ){ |
736 |
|
737 |
if( readLine[0] != '!' ){ |
738 |
|
739 |
the_token = strtok( readLine, " \n\t,;" ); |
740 |
if( the_token != NULL ){ |
741 |
|
742 |
strcpy( currentTorsionType->nameA, the_token ); |
743 |
|
744 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
745 |
fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum ); |
746 |
exit(8); |
747 |
} |
748 |
|
749 |
strcpy( currentTorsionType->nameB, the_token ); |
750 |
|
751 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
752 |
fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum ); |
753 |
exit(8); |
754 |
} |
755 |
|
756 |
strcpy( currentTorsionType->nameC, the_token ); |
757 |
|
758 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
759 |
fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum ); |
760 |
exit(8); |
761 |
} |
762 |
|
763 |
strcpy( currentTorsionType->nameD, the_token ); |
764 |
|
765 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
766 |
fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum ); |
767 |
exit(8); |
768 |
} |
769 |
|
770 |
strcpy( currentTorsionType->type, the_token ); |
771 |
|
772 |
if( !strcmp( currentTorsionType->type, "cubic" ) ){ |
773 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
774 |
fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum ); |
775 |
exit(8); |
776 |
} |
777 |
|
778 |
sscanf( the_token, "%lf", ¤tTorsionType->k1 ); |
779 |
|
780 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
781 |
fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum ); |
782 |
exit(8); |
783 |
} |
784 |
|
785 |
sscanf( the_token, "%lf", ¤tTorsionType->k2 ); |
786 |
|
787 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
788 |
fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum ); |
789 |
exit(8); |
790 |
} |
791 |
|
792 |
sscanf( the_token, "%lf", ¤tTorsionType->k3 ); |
793 |
|
794 |
if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ |
795 |
fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum ); |
796 |
exit(8); |
797 |
} |
798 |
|
799 |
sscanf( the_token, "%lf", ¤tTorsionType->k4 ); |
800 |
} |
801 |
|
802 |
else{ |
803 |
fprintf(stderr, |
804 |
"Unknown TraPPE torsion type \"%s\" at line %d\n", |
805 |
currentTorsionType->type, |
806 |
lineNum ); |
807 |
exit(8); |
808 |
} |
809 |
} |
810 |
} |
811 |
|
812 |
tempTorsionType = new LinkedType; |
813 |
currentTorsionType->next = tempTorsionType; |
814 |
currentTorsionType = tempTorsionType; |
815 |
|
816 |
eof_test = fgets( readLine, sizeof(readLine), frcFile ); |
817 |
lineNum++; |
818 |
} |
819 |
|
820 |
|
821 |
// initialize the Torsions |
822 |
|
823 |
for( i=0; i<nTorsions; i++ ){ |
824 |
|
825 |
a = the_torsions[i].a; |
826 |
b = the_torsions[i].b; |
827 |
c = the_torsions[i].c; |
828 |
d = the_torsions[i].d; |
829 |
|
830 |
atomA = the_atoms[a]->getType(); |
831 |
atomB = the_atoms[b]->getType(); |
832 |
atomC = the_atoms[c]->getType(); |
833 |
atomD = the_atoms[d]->getType(); |
834 |
currentTorsionType = headTorsionType->find( atomA, atomB, atomC, atomD ); |
835 |
if( currentTorsionType == NULL ){ |
836 |
fprintf( stderr, "TorsionType error, %s - %s - %s - %s not found" |
837 |
" in force file.\n", |
838 |
atomA, atomB, atomC, atomD ); |
839 |
exit(8); |
840 |
} |
841 |
|
842 |
if( !strcmp( currentTorsionType->type, "cubic" ) ){ |
843 |
index = i + entry_plug->n_bonds + entry_plug->n_bends; |
844 |
|
845 |
cTors = new CubicTorsion( *the_atoms[a], *the_atoms[b], |
846 |
*the_atoms[c], *the_atoms[d] ); |
847 |
cTors->setConstants( currentTorsionType->k1, currentTorsionType->k2, |
848 |
currentTorsionType->k3, currentTorsionType->k4 ); |
849 |
the_sris[index] = cTors; |
850 |
} |
851 |
} |
852 |
|
853 |
|
854 |
// clean up the memory |
855 |
|
856 |
delete headTorsionType; |
857 |
} |