1
2 #include "msi2lmp.h"
3 #include "Forcefield.h"
4
5 #include <string.h>
6 #include <stdlib.h>
7 #include <math.h>
8
9 static int find_improper_body_data(char [][5],struct FrcFieldItem,int *);
10 static void rearrange_improper(int,int);
11 static int find_trigonal_body_data(char [][5],struct FrcFieldItem);
12 static int find_angleangle_data(char [][5],struct FrcFieldItem,int[]);
13 static int find_match(int, char [][5],struct FrcFieldItem,int *);
14 static int match_types(int,int,char [][5],char [][5],int *);
15 static double get_r0(int,int);
16 static double get_t0(int,int,int);
17 static int quo_cp();
18 static void get_equivs(int,char [][5],char[][5]);
19 static int find_equiv_type(char[]);
20
21 /**********************************************************************/
22 /* */
23 /* GetParameters is a long routine for searching the forcefield */
24 /* parameters (read in by ReadFrcFile) for parameters corresponding */
25 /* to the different internal coordinate types derived by MakeLists */
26 /* */
27 /**********************************************************************/
28
GetParameters()29 void GetParameters()
30 {
31 int i,j,k,backwards,cp_type,rearrange;
32 int kloc[3],multiplicity;
33 char potential_types[4][5];
34 char equiv_types[4][5];
35 double rab,rbc,rcd,tabc,tbcd,tabd,tcbd;
36
37 if (pflag > 1) fprintf(stderr," Trying Atom Equivalences if needed\n");
38
39 /**********************************************************************/
40 /* */
41 /* Find masses of atom types */
42 /* */
43 /**********************************************************************/
44
45 for (i=0; i < no_atom_types; i++) {
46 backwards = -1;
47 strncpy(potential_types[0],atomtypes[i].potential,5);
48 k = find_match(1,potential_types,ff_atomtypes,&backwards);
49 if (k < 0) {
50 printf(" Unable to find mass for %s\n",atomtypes[i].potential);
51 condexit(10);
52 } else {
53 atomtypes[i].mass = ff_atomtypes.data[k].ff_param[0];
54 }
55 }
56
57 /**********************************************************************/
58 /* */
59 /* Find VDW parameters for atom types */
60 /* */
61 /**********************************************************************/
62
63 for (i=0; i < no_atom_types; i++) {
64 backwards = 0;
65 for (j=0; j < 2; j++) atomtypes[i].params[j] = 0.0;
66 strncpy(potential_types[0],atomtypes[i].potential,5);
67 k = find_match(1,potential_types,ff_vdw,&backwards);
68 if (k < 0) {
69 get_equivs(1,potential_types,equiv_types);
70
71 if (pflag > 2) printf(" Using equivalences for VDW %s -> %s\n",
72 potential_types[0],equiv_types[0]);
73
74 k = find_match(1,equiv_types,ff_vdw,&backwards);
75 }
76 if (k < 0) {
77 printf(" Unable to find vdw data for %s\n",atomtypes[i].potential);
78 condexit(11);
79 } else {
80 if (ljtypeflag == 0) {
81 if((ff_vdw.data[k].ff_param[0] != 0.0 ) &&
82 (ff_vdw.data[k].ff_param[1] != 0.0)) {
83 atomtypes[i].params[0] =
84 (ff_vdw.data[k].ff_param[1]*
85 ff_vdw.data[k].ff_param[1])/(4.0*ff_vdw.data[k].ff_param[0]);
86 atomtypes[i].params[1] = pow((ff_vdw.data[k].ff_param[0]/
87 ff_vdw.data[k].ff_param[1]),
88 (1.0/6.0));
89 }
90 } else if (ljtypeflag == 1) {
91 atomtypes[i].params[0] = ff_vdw.data[k].ff_param[1];
92 atomtypes[i].params[1] = ff_vdw.data[k].ff_param[0];
93 } else {
94 printf(" Unknown LJ parameter type %d\n",ljtypeflag);
95 exit(111);
96 }
97 }
98 }
99
100 if (pflag > 2) {
101 printf("\n Atom Types, Masses and VDW Parameters\n");
102 for (i=0; i < no_atom_types; i++) {
103 printf(" %3s %8.4f %8.4f %8.4f\n",
104 atomtypes[i].potential,atomtypes[i].mass, atomtypes[i].params[0],atomtypes[i].params[1]);
105 }
106 }
107
108 /**********************************************************************/
109 /* */
110 /* Find parameters for bond types */
111 /* */
112 /**********************************************************************/
113
114 for (i=0; i < no_bond_types; i++) {
115 backwards = 0;
116 for (j=0; j < 4; j++) bondtypes[i].params[j] = 0.0;
117 for (j=0; j < 2; j++)
118 strncpy(potential_types[j],
119 atomtypes[bondtypes[i].types[j]].potential,5);
120 k = find_match(2,potential_types,ff_bond,&backwards);
121 if (k < 0) {
122 get_equivs(2,potential_types,equiv_types);
123
124 if (pflag > 2) {
125 printf(" Using equivalences for bond %s %s -> %s %s\n",
126 potential_types[0],potential_types[1],
127 equiv_types[0],equiv_types[1]);
128 }
129 k = find_match(2,equiv_types,ff_bond,&backwards);
130 }
131 if (k < 0) {
132 printf(" Unable to find bond data for %s %s\n",
133 potential_types[0],potential_types[1]);
134 condexit(12);
135 } else {
136 if (forcefield & (FF_TYPE_CLASS1|FF_TYPE_OPLSAA)) {
137 bondtypes[i].params[0] = ff_bond.data[k].ff_param[1];
138 bondtypes[i].params[1] = ff_bond.data[k].ff_param[0];
139 }
140
141 if (forcefield & FF_TYPE_CLASS2) {
142 for (j=0; j < 4; j++)
143 bondtypes[i].params[j] = ff_bond.data[k].ff_param[j];
144 }
145 }
146 }
147
148 if (pflag > 2) {
149 printf("\n Bond Types and Parameters\n");
150 for (i=0; i < no_bond_types; i++) {
151 for (j=0; j < 2; j++)
152 printf(" %-3s",atomtypes[bondtypes[i].types[j]].potential);
153 for (j=0; j < 4; j++)
154 printf(" %8.4f",bondtypes[i].params[j]);
155 printf("\n");
156 }
157 }
158
159
160 /**********************************************************************/
161 /* */
162 /* Find parameters for angle types including bondbond, */
163 /* and bondangle parameters if Class II */
164 /* */
165 /* Each of the cross terms are searched separately even though */
166 /* they share a given angle type. This allows parameters to be */
167 /* in different order in the forcefield for each cross term or */
168 /* maybe not even there. */
169 /* */
170 /**********************************************************************/
171 for (i=0; i < no_angle_types; i++) {
172 backwards = 0;
173 for (j=0; j < 4; j++) angletypes[i].params[j] = 0.0;
174 for (j=0; j < 3; j++)
175 strncpy(potential_types[j],atomtypes[angletypes[i].types[j]].potential,5);
176 k = find_match(3,potential_types,ff_ang,&backwards);
177 if (k < 0) {
178 get_equivs(3,potential_types,equiv_types);
179 if (pflag > 2) {
180 printf(" Using equivalences for angle %s %s %s -> %s %s %s\n",
181 potential_types[0],potential_types[1],
182 potential_types[2],
183 equiv_types[0],equiv_types[1],
184 equiv_types[2]);
185 }
186 k = find_match(3,equiv_types,ff_ang,&backwards);
187 }
188 if (k < 0) {
189 printf(" Unable to find angle data for %s %s %s\n",
190 potential_types[0],potential_types[1],potential_types[2]);
191 condexit(13);
192 } else {
193 if (forcefield & (FF_TYPE_CLASS1|FF_TYPE_OPLSAA)) {
194 angletypes[i].params[0] = ff_ang.data[k].ff_param[1];
195 angletypes[i].params[1] = ff_ang.data[k].ff_param[0];
196 }
197
198 if (forcefield & FF_TYPE_CLASS2) {
199 for (j=0; j < 4; j++)
200 angletypes[i].params[j] = ff_ang.data[k].ff_param[j];
201 }
202 }
203 if (forcefield & FF_TYPE_CLASS2) {
204 get_equivs(3,potential_types,equiv_types);
205 if (pflag > 2) {
206 printf(" Using equivalences for 3 body cross terms %s %s %s -> %s %s %s\n",
207 potential_types[0],potential_types[1],potential_types[2],
208 equiv_types[0],equiv_types[1],equiv_types[2]);
209 }
210 for (j=0; j < 3; j++) angletypes[i].bondbond_cross_term[j] = 0.0;
211 for (j=0; j < 4; j++) angletypes[i].bondangle_cross_term[j] = 0.0;
212
213 rab = get_r0(angletypes[i].types[0],angletypes[i].types[1]);
214 rbc = get_r0(angletypes[i].types[1],angletypes[i].types[2]);
215
216 angletypes[i].bondbond_cross_term[1] = rab;
217 angletypes[i].bondbond_cross_term[2] = rbc;
218 angletypes[i].bondangle_cross_term[2] = rab;
219 angletypes[i].bondangle_cross_term[3] = rbc;
220
221 k = find_match(3,potential_types,ff_bonbon,&backwards);
222 if (k < 0) {
223 k = find_match(3,equiv_types,ff_bonbon,&backwards);
224 }
225 if (k < 0) {
226 printf(" Unable to find bondbond data for %s %s %s\n",
227 potential_types[0],potential_types[1],potential_types[2]);
228 condexit(14);
229 } else {
230 angletypes[i].bondbond_cross_term[0] = ff_bonbon.data[k].ff_param[0];
231 }
232 k = find_match(3,potential_types,ff_bonang,&backwards);
233 if (k < 0) {
234 k = find_match(3,equiv_types,ff_bonang,&backwards);
235 }
236 if (k < 0) {
237 printf(" Unable to find bondangle data for %s %s %s\n",
238 potential_types[0],potential_types[1],potential_types[2]);
239 condexit(15);
240 } else {
241 if (backwards) {
242 angletypes[i].bondangle_cross_term[0] = ff_bonang.data[k].ff_param[1];
243 angletypes[i].bondangle_cross_term[1] = ff_bonang.data[k].ff_param[0];
244 } else {
245 angletypes[i].bondangle_cross_term[0] = ff_bonang.data[k].ff_param[0];
246 angletypes[i].bondangle_cross_term[1] = ff_bonang.data[k].ff_param[1];
247 }
248 }
249 }
250 }
251
252 if (pflag > 2) {
253 printf("\n Angle Types and Parameters\n");
254 for (i=0; i < no_angle_types; i++) {
255 for (j=0; j < 3; j++)
256 printf(" %-3s", atomtypes[angletypes[i].types[j]].potential);
257 for (j=0; j < 4; j++) printf(" %8.4f",angletypes[i].params[j]);
258 printf("\n");
259 }
260
261 if (forcefield & FF_TYPE_CLASS2) {
262 printf("\n BondBond Types and Parameters\n");
263 for (i=0; i < no_angle_types; i++) {
264 for (j=0; j < 3; j++)
265 printf(" %-3s",atomtypes[angletypes[i].types[j]].potential);
266 for (j=0; j < 3; j++)
267 printf(" %8.4f",angletypes[i].bondbond_cross_term[j]);
268 printf("\n");
269 }
270 printf("\n BondAngle Types and Parameters\n");
271 for (i=0; i < no_angle_types; i++) {
272 for (j=0; j < 3; j++)
273 printf(" %-3s",atomtypes[angletypes[i].types[j]].potential);
274 for (j=0; j < 4; j++)
275 printf(" %8.4f",angletypes[i].bondangle_cross_term[j]);
276 printf("\n");
277 }
278 }
279 }
280
281 /**********************************************************************/
282 /* */
283 /* Find parameters for dihedral types including endbonddihedral, */
284 /* midbonddihedral, angledihedral, angleangledihedral and */
285 /* bondbond13 parameters if Class II */
286 /* */
287 /* Each of the cross terms are searched separately even though */
288 /* they share a given dihedral type. This allows parameters to be */
289 /* in different order in the forcefield for each cross term or */
290 /* maybe not even there. */
291 /* */
292 /**********************************************************************/
293
294 for (i=0; i < no_dihedral_types; i++) {
295 for (j=0; j < 6; j++)
296 dihedraltypes[i].params[j] = 0.0;
297 for (j=0; j < 4; j++)
298 strncpy(potential_types[j],
299 atomtypes[dihedraltypes[i].types[j]].potential,5);
300 backwards = 0;
301 k = find_match(4,potential_types,ff_tor,&backwards);
302
303 if (k < 0) {
304 get_equivs(4,potential_types,equiv_types);
305
306 if (pflag > 2) {
307 printf(" Using equivalences for dihedral %s %s %s %s -> %s %s %s %s\n",
308 potential_types[0],potential_types[1],
309 potential_types[2],potential_types[3],
310 equiv_types[0],equiv_types[1],
311 equiv_types[2],equiv_types[3]);
312 }
313 k = find_match(4,equiv_types,ff_tor,&backwards);
314 }
315 if (k < 0) {
316 printf(" Unable to find torsion data for %s %s %s %s\n",
317 potential_types[0],
318 potential_types[1],
319 potential_types[2],
320 potential_types[3]);
321 condexit(16);
322 } else {
323 if (forcefield & FF_TYPE_CLASS1) {
324 multiplicity = 1;
325 if (ff_tor.data[k].ff_types[0][0] == '*')
326 multiplicity =
327 atomtypes[dihedraltypes[i].types[1]].no_connect-1;
328 if (ff_tor.data[k].ff_types[3][0] == '*')
329 multiplicity *=
330 atomtypes[dihedraltypes[i].types[2]].no_connect-1;
331
332 dihedraltypes[i].params[0] = ff_tor.data[k].ff_param[0]/(double) multiplicity;
333 if (ff_tor.data[k].ff_param[2] == 0.0)
334 dihedraltypes[i].params[1] = 1.0;
335 else if (ff_tor.data[k].ff_param[2] == 180.0)
336 dihedraltypes[i].params[1] = -1.0;
337 else {
338 printf(" Non planar phi0 for %s %s %s %s\n",
339 potential_types[0],potential_types[1],
340 potential_types[2],potential_types[3]);
341 dihedraltypes[i].params[1] = 0.0;
342 }
343 dihedraltypes[i].params[2] = ff_tor.data[k].ff_param[1];
344 }
345 if (forcefield & FF_TYPE_OPLSAA) {
346 for (j=0; j < 4; j++)
347 dihedraltypes[i].params[j] = ff_tor.data[k].ff_param[j];
348 }
349 if (forcefield & FF_TYPE_CLASS2) {
350 for (j=0; j < 6; j++)
351 dihedraltypes[i].params[j] = ff_tor.data[k].ff_param[j];
352 }
353 }
354
355 if (forcefield & FF_TYPE_CLASS2) {
356 get_equivs(4,potential_types,equiv_types);
357 if (pflag > 2) {
358 printf(" Using equivalences for linear 4 body cross terms %s %s %s %s -> %s %s %s %s\n",
359 potential_types[0],potential_types[1],
360 potential_types[2],potential_types[3],
361 equiv_types[0],equiv_types[1],
362 equiv_types[2],equiv_types[3]);
363 }
364
365 for (j=0; j < 8; j++)
366 dihedraltypes[i].endbonddihedral_cross_term[j] = 0.0;
367 for (j=0; j < 4; j++)
368 dihedraltypes[i].midbonddihedral_cross_term[j] = 0.0;
369 for (j=0; j < 8; j++)
370 dihedraltypes[i].angledihedral_cross_term[j] = 0.0;
371 for (j=0; j < 3; j++)
372 dihedraltypes[i].angleangledihedral_cross_term[j] = 0.0;
373 for (j=0; j < 3; j++)
374 dihedraltypes[i].bond13_cross_term[j] = 0.0;
375
376 rab = get_r0(dihedraltypes[i].types[0],dihedraltypes[i].types[1]);
377 rbc = get_r0(dihedraltypes[i].types[1],dihedraltypes[i].types[2]);
378 rcd = get_r0(dihedraltypes[i].types[2],dihedraltypes[i].types[3]);
379 tabc = get_t0(dihedraltypes[i].types[0],
380 dihedraltypes[i].types[1],
381 dihedraltypes[i].types[2]);
382
383 tbcd = get_t0(dihedraltypes[i].types[1],
384 dihedraltypes[i].types[2],
385 dihedraltypes[i].types[3]);
386
387 dihedraltypes[i].endbonddihedral_cross_term[6] = rab;
388 dihedraltypes[i].endbonddihedral_cross_term[7] = rcd;
389 dihedraltypes[i].midbonddihedral_cross_term[3] = rbc;
390 dihedraltypes[i].angledihedral_cross_term[6] = tabc;
391 dihedraltypes[i].angledihedral_cross_term[7] = tbcd;
392 dihedraltypes[i].angleangledihedral_cross_term[1] = tabc;
393 dihedraltypes[i].angleangledihedral_cross_term[2] = tbcd;
394 dihedraltypes[i].bond13_cross_term[1] = rab;
395 dihedraltypes[i].bond13_cross_term[2] = rcd;
396
397 backwards = 0;
398 k = find_match(4,potential_types,ff_endbontor,&backwards);
399 if (k < 0) {
400 k = find_match(4,equiv_types,ff_endbontor,&backwards);
401 }
402 if (k < 0) {
403 printf(" Unable to find endbonddihedral data for %s %s %s %s\n",
404 potential_types[0],potential_types[1],
405 potential_types[2],potential_types[3]);
406 condexit(17);
407 } else {
408 if (backwards) {
409 dihedraltypes[i].endbonddihedral_cross_term[0] =
410 ff_endbontor.data[k].ff_param[3];
411 dihedraltypes[i].endbonddihedral_cross_term[1] =
412 ff_endbontor.data[k].ff_param[4];
413 dihedraltypes[i].endbonddihedral_cross_term[2] =
414 ff_endbontor.data[k].ff_param[5];
415 dihedraltypes[i].endbonddihedral_cross_term[3] =
416 ff_endbontor.data[k].ff_param[0];
417 dihedraltypes[i].endbonddihedral_cross_term[4] =
418 ff_endbontor.data[k].ff_param[1];
419 dihedraltypes[i].endbonddihedral_cross_term[5] =
420 ff_endbontor.data[k].ff_param[2];
421 } else {
422 dihedraltypes[i].endbonddihedral_cross_term[0] =
423 ff_endbontor.data[k].ff_param[0];
424 dihedraltypes[i].endbonddihedral_cross_term[1] =
425 ff_endbontor.data[k].ff_param[1];
426 dihedraltypes[i].endbonddihedral_cross_term[2] =
427 ff_endbontor.data[k].ff_param[2];
428 dihedraltypes[i].endbonddihedral_cross_term[3] =
429 ff_endbontor.data[k].ff_param[3];
430 dihedraltypes[i].endbonddihedral_cross_term[4] =
431 ff_endbontor.data[k].ff_param[4];
432 dihedraltypes[i].endbonddihedral_cross_term[5] =
433 ff_endbontor.data[k].ff_param[5];
434 }
435 }
436 backwards = 0;
437 k = find_match(4,potential_types,ff_midbontor,&backwards);
438 if (k < 0) {
439 k = find_match(4,equiv_types,ff_midbontor,&backwards);
440 }
441 if (k < 0) {
442 printf(" Unable to find midbonddihedral data for %s %s %s %s\n",
443 potential_types[0],potential_types[1],
444 potential_types[2],potential_types[3]);
445 condexit(18);
446 } else {
447 dihedraltypes[i].midbonddihedral_cross_term[0] =
448 ff_midbontor.data[k].ff_param[0];
449 dihedraltypes[i].midbonddihedral_cross_term[1] =
450 ff_midbontor.data[k].ff_param[1];
451 dihedraltypes[i].midbonddihedral_cross_term[2] =
452 ff_midbontor.data[k].ff_param[2];
453 }
454
455 backwards = 0;
456 k = find_match(4,potential_types,ff_angtor,&backwards);
457 if (k < 0) {
458 k = find_match(4,equiv_types,ff_angtor,&backwards);
459 }
460 if (k < 0) {
461 printf(" Unable to find angledihedral data for %s %s %s %s\n",
462 potential_types[0],potential_types[1],
463 potential_types[2],potential_types[3]);
464 condexit(19);
465 } else {
466 if (backwards) {
467 dihedraltypes[i].angledihedral_cross_term[0] =
468 ff_angtor.data[k].ff_param[3];
469 dihedraltypes[i].angledihedral_cross_term[1] =
470 ff_angtor.data[k].ff_param[4];
471 dihedraltypes[i].angledihedral_cross_term[2] =
472 ff_angtor.data[k].ff_param[5];
473 dihedraltypes[i].angledihedral_cross_term[3] =
474 ff_angtor.data[k].ff_param[0];
475 dihedraltypes[i].angledihedral_cross_term[4] =
476 ff_angtor.data[k].ff_param[1];
477 dihedraltypes[i].angledihedral_cross_term[5] =
478 ff_angtor.data[k].ff_param[2];
479 } else {
480 dihedraltypes[i].angledihedral_cross_term[0] =
481 ff_angtor.data[k].ff_param[0];
482 dihedraltypes[i].angledihedral_cross_term[1] =
483 ff_angtor.data[k].ff_param[1];
484 dihedraltypes[i].angledihedral_cross_term[2] =
485 ff_angtor.data[k].ff_param[2];
486 dihedraltypes[i].angledihedral_cross_term[3] =
487 ff_angtor.data[k].ff_param[3];
488 dihedraltypes[i].angledihedral_cross_term[4] =
489 ff_angtor.data[k].ff_param[4];
490 dihedraltypes[i].angledihedral_cross_term[5] =
491 ff_angtor.data[k].ff_param[5];
492 }
493 }
494 backwards = 0;
495 k = find_match(4,potential_types,ff_angangtor,&backwards);
496 if (k < 0) {
497 k = find_match(4,equiv_types,ff_angangtor,&backwards);
498 }
499 if (k < 0) {
500 printf(" Unable to find angleangledihedral data for %s %s %s %s\n",
501 potential_types[0],potential_types[1],
502 potential_types[2],potential_types[3]);
503 condexit(20);
504 } else {
505 dihedraltypes[i].angleangledihedral_cross_term[0] =
506 ff_angangtor.data[k].ff_param[0];
507 }
508 cp_type = quo_cp();
509 if ((cp_type >= 0) &&
510 ((dihedraltypes[i].types[0] == cp_type) ||
511 (dihedraltypes[i].types[1] == cp_type) ||
512 (dihedraltypes[i].types[2] == cp_type) ||
513 (dihedraltypes[i].types[3] == cp_type) )) {
514 backwards = 0;
515 k = find_match(4,potential_types,ff_bonbon13,&backwards);
516 if (k < 0) {
517 k = find_match(4,equiv_types,ff_bonbon13,&backwards);
518 }
519 if (k < 0) {
520 printf(" Unable to find bond13 data for %s %s %s %s\n",
521 potential_types[0],potential_types[1],
522 potential_types[2],potential_types[3]);
523 condexit(21);
524 } else {
525 dihedraltypes[i].bond13_cross_term[0] =
526 ff_bonbon13.data[k].ff_param[0];
527 }
528 }
529 }
530 }
531
532 if (pflag > 2) {
533 printf("\n Dihedral Types and Parameters\n");
534 for (i=0; i < no_dihedral_types; i++) {
535 for (j=0; j < 4; j++)
536 printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential);
537 for (j=0; j < 6; j++)
538 printf(" %8.4f",dihedraltypes[i].params[j]);
539 printf("\n");
540 }
541
542 if (forcefield & FF_TYPE_CLASS2) {
543
544 printf("\n EndBondDihedral Types and Parameters\n");
545 for (i=0; i < no_dihedral_types; i++) {
546 for (j=0; j < 4; j++)
547 printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential);
548 for (j=0; j < 8; j++)
549 printf(" %8.4f",dihedraltypes[i].endbonddihedral_cross_term[j]);
550 printf("\n");
551 }
552 printf("\n MidBondDihedral Types and Parameters\n");
553 for (i=0; i < no_dihedral_types; i++) {
554 for (j=0; j < 4; j++)
555 printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential);
556 for (j=0; j < 4; j++)
557 printf(" %8.4f",dihedraltypes[i].midbonddihedral_cross_term[j]);
558 printf("\n");
559 }
560
561 printf("\n AngleDihedral Types and Parameters\n");
562 for (i=0; i < no_dihedral_types; i++) {
563 for (j=0; j < 4; j++)
564 printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential);
565 for (j=0; j < 8; j++)
566 printf(" %8.4f",dihedraltypes[i].angledihedral_cross_term[j]);
567 printf("\n");
568 }
569
570 printf("\n AngleAngleDihedral Types and Parameters\n");
571 for (i=0; i < no_dihedral_types; i++) {
572 for (j=0; j < 4; j++)
573 printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential);
574 for (j=0; j < 3; j++)
575 printf(" %8.4f",dihedraltypes[i].angleangledihedral_cross_term[j]);
576 printf("\n");
577 }
578
579 printf("\n Bond13 Types and Parameters\n");
580
581 for (i=0; i < no_dihedral_types; i++) {
582 for (j=0; j < 4; j++)
583 printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential);
584 for (j=0; j < 3; j++)
585 printf(" %8.4f",dihedraltypes[i].bond13_cross_term[j]);
586 printf("\n");
587 }
588 }
589 }
590
591
592 /**********************************************************************/
593 /* */
594 /* Find parameters for oop types */
595 /* */
596 /* This is the most complicated of all the types because the */
597 /* the class I oop is actually an improper torsion and does */
598 /* not have the permutation symmetry of a well defined oop */
599 /* The net result is that if one does not find the current */
600 /* atom type ordering in the forcefield file then one must try each */
601 /* of the next permutations (6 in total) and when a match is found */
602 /* the program must go back and rearrange the oop type AND the atom */
603 /* ordering in the oop lists for those with the current type */
604 /* */
605 /* The Class II oop types are easier but also tedious since the */
606 /* program has to try all permutations of the a c and d atom */
607 /* types to find a match. A special routine is used to do this. */
608 /* */
609 /* Fortunately, there are typically few oop types */
610 /* */
611 /**********************************************************************/
612
613 if (forcefield & FF_TYPE_CLASS1) {
614 for (i=0; i < no_oop_types; i++) {
615 for (j=0; j < 3; j++) ooptypes[i].params[j] = 0.0;
616 for (j=0; j < 4; j++)
617 strncpy(potential_types[j],
618 atomtypes[ooptypes[i].types[j]].potential,5);
619
620 k = find_improper_body_data(potential_types,ff_oop,&rearrange);
621 if (k < 0) {
622 get_equivs(5,potential_types,equiv_types);
623
624 if (pflag > 2) {
625 printf(" Using equivalences for oop %s %s %s %s -> %s %s %s %s\n",
626 potential_types[0],potential_types[1],
627 potential_types[2],potential_types[3],
628 equiv_types[0],equiv_types[1],
629 equiv_types[2],equiv_types[3]);
630 }
631 k = find_improper_body_data(equiv_types,ff_oop,&rearrange);
632 }
633 if (k < 0) {
634 printf(" Unable to find oop data for %s %s %s %s\n",
635 potential_types[0],
636 potential_types[1],potential_types[2],potential_types[3]);
637 condexit(22);
638 } else {
639 ooptypes[i].params[0] = ff_oop.data[k].ff_param[0];
640 if (ff_oop.data[k].ff_param[2] == 0.0)
641 ooptypes[i].params[1] = 1.0;
642 else if (ff_oop.data[k].ff_param[2] == 180.0)
643 ooptypes[i].params[1] = -1.0;
644 else {
645 printf(" Non planar phi0 for %s %s %s %s\n",
646 potential_types[0],potential_types[1],
647 potential_types[2],potential_types[3]);
648 ooptypes[i].params[1] = 0.0;
649 }
650 ooptypes[i].params[2] = ff_oop.data[k].ff_param[1];
651 if (rearrange > 0) rearrange_improper(i,rearrange);
652 }
653 }
654 }
655
656 if (forcefield & FF_TYPE_CLASS2) {
657 for (i=0; i < no_oop_types; i++) {
658 for (j=0; j < 3; j++)
659 ooptypes[i].params[j] = 0.0;
660 for (j=0; j < 4; j++)
661 strncpy(potential_types[j],
662 atomtypes[ooptypes[i].types[j]].potential,5);
663 k = find_trigonal_body_data(potential_types,ff_oop);
664 if (k < 0) {
665 get_equivs(5,potential_types,equiv_types);
666 if (pflag > 2) {
667 printf(" Using equivalences for oop %s %s %s %s -> %s %s %s %s\n",
668 potential_types[0],potential_types[1],
669 potential_types[2],potential_types[3],
670 equiv_types[0],equiv_types[1],
671 equiv_types[2],equiv_types[3]);
672 }
673 k = find_trigonal_body_data(equiv_types,ff_oop);
674 }
675 if (k < 0) {
676 printf(" Unable to find oop data for %s %s %s %s\n",
677 potential_types[0],
678 potential_types[1],potential_types[2],potential_types[3]);
679 condexit(23);
680 } else {
681 for (j=0; j < 2; j++)
682 ooptypes[i].params[j] = ff_oop.data[k].ff_param[j];
683 }
684 }
685 }
686
687 if (pflag > 2) {
688 printf("\n OOP Types and Parameters\n");
689 for (i=0; i < no_oop_types; i++) {
690 for (j=0; j < 4; j++)
691 printf(" %-3s",atomtypes[ooptypes[i].types[j]].potential);
692 for (j=0; j < 3; j++)
693 printf(" %8.4f",ooptypes[i].params[j]);
694 printf("\n");
695 }
696 }
697
698
699 /**********************************************************************/
700 /* */
701 /* Find parameters for angleangle types (Class II only) */
702 /* */
703 /* This is somewhat complicated in that one set of four types */
704 /* a b c d has three angleangle combinations so for each type */
705 /* the program needs to find three sets of parameters by */
706 /* progressively looking for data for different permutations of */
707 /* a c and d */
708 /* */
709 /**********************************************************************/
710
711 if (forcefield & FF_TYPE_CLASS2) {
712
713 for (i=0; i < no_oop_types; i++) {
714
715 for (j=0; j < 6; j++) ooptypes[i].angleangle_params[j] = 0.0;
716
717 for (j=0; j < 4; j++)
718 strncpy(potential_types[j],
719 atomtypes[ooptypes[i].types[j]].potential,5);
720
721
722 tabc = get_t0(ooptypes[i].types[0],
723 ooptypes[i].types[1],
724 ooptypes[i].types[2]);
725
726 tabd = get_t0(ooptypes[i].types[0],
727 ooptypes[i].types[1],
728 ooptypes[i].types[3]);
729 tcbd = get_t0(ooptypes[i].types[2],
730 ooptypes[i].types[1],
731
732 ooptypes[i].types[3]);
733
734 ooptypes[i].angleangle_params[3] = tabc;
735 ooptypes[i].angleangle_params[4] = tcbd;
736 ooptypes[i].angleangle_params[5] = tabd;
737
738 k = find_angleangle_data(potential_types,ff_angang,kloc);
739 if (k < 0) {
740 get_equivs(5,potential_types,equiv_types);
741 if (pflag > 2) {
742 printf(" Using equivalences for angleangle %s %s %s %s -> %s %s %s %s\n",
743 potential_types[0],potential_types[1],
744 potential_types[2],potential_types[3],
745 equiv_types[0],equiv_types[1],
746 equiv_types[2],equiv_types[3]);
747 k = find_angleangle_data(equiv_types,ff_angang,kloc);
748 }
749 }
750 if (k < 0) {
751 printf(" Unable to find angleangle data for %s %s %s %s\n",
752 potential_types[0],
753 potential_types[1],potential_types[2],potential_types[3]);
754 condexit(24);
755 } else {
756 for (j=0; j < 3; j++) {
757 if (kloc[j] > -1)
758 ooptypes[i].angleangle_params[j] = ff_angang.data[kloc[j]].ff_param[0];
759 }
760 }
761 }
762
763 for (i=0; i < no_angleangle_types; i++) {
764 for (j=0; j < 6; j++) angleangletypes[i].params[j] = 0.0;
765 for (j=0; j < 4; j++)
766 strncpy(potential_types[j],
767 atomtypes[angleangletypes[i].types[j]].potential,5);
768
769 tabc = get_t0(angleangletypes[i].types[0],
770 angleangletypes[i].types[1],
771 angleangletypes[i].types[2]);
772 tabd = get_t0(angleangletypes[i].types[0],
773 angleangletypes[i].types[1],
774 angleangletypes[i].types[3]);
775 tcbd = get_t0(angleangletypes[i].types[2],
776 angleangletypes[i].types[1],
777 angleangletypes[i].types[3]);
778
779 angleangletypes[i].params[3] = tabc;
780 angleangletypes[i].params[4] = tcbd;
781 angleangletypes[i].params[5] = tabd;
782
783 k = find_angleangle_data(potential_types,ff_angang,kloc);
784 if (k < 0) {
785 get_equivs(5,potential_types,equiv_types);
786 if (pflag > 2) {
787 printf("Using equivalences for angleangle %s %s %s %s -> %s %s %s %s\n",
788 potential_types[0],potential_types[1],
789 potential_types[2],potential_types[3],
790 equiv_types[0],equiv_types[1],
791 equiv_types[2],equiv_types[3]);
792 }
793 k = find_angleangle_data(equiv_types,ff_angang,kloc);
794 }
795 if (k < 0) {
796 printf(" Unable to find angleangle data for %s %s %s %s\n",
797 potential_types[0],
798 potential_types[1],potential_types[2],potential_types[3]);
799 condexit(25);
800 } else {
801 for (j=0; j < 3; j++) {
802 if (kloc[j] > -1)
803 angleangletypes[i].params[j] =
804 ff_angang.data[kloc[j]].ff_param[0];
805 }
806 }
807 }
808 if (pflag > 2) {
809 printf("\n AngleAngle Types and Parameters\n");
810 for (i=0; i < no_oop_types; i++) {
811 for (j=0; j < 4; j++)
812 printf(" %-3s",atomtypes[ooptypes[i].types[j]].potential);
813 for (j=0; j < 6; j++)
814 printf(" %8.4f",ooptypes[i].angleangle_params[j]);
815 printf("\n");
816 }
817 for (i=0; i < no_angleangle_types; i++) {
818 for (j=0; j < 4; j++)
819 printf(" %-3s",atomtypes[angleangletypes[i].types[j]].potential);
820 for (j=0; j < 6; j++) printf(" %8.4f",angleangletypes[i].params[j]);
821 printf("\n");
822 }
823 }
824 }
825 }
826
find_improper_body_data(char types1[][5],struct FrcFieldItem item,int * rearrange_ptr)827 int find_improper_body_data(char types1[][5],struct FrcFieldItem item,
828 int *rearrange_ptr)
829 {
830 int k,backwards;
831 char mirror_types[4][5];
832
833 backwards = 0;
834
835 /* a b c d */
836
837 *rearrange_ptr = 0;
838 k = find_match(4,types1,item,&backwards);
839 if (k >= 0) return k;
840
841 /* a b d c */
842
843 *rearrange_ptr = 1;
844 strncpy(mirror_types[0],types1[0],5);
845 strncpy(mirror_types[1],types1[1],5);
846 strncpy(mirror_types[2],types1[3],5);
847 strncpy(mirror_types[3],types1[2],5);
848 k = find_match(4,mirror_types,item,&backwards);
849 if (k >= 0) return k;
850
851 /* d b a c */
852
853 *rearrange_ptr = 2;
854 strncpy(mirror_types[0],types1[3],5);
855 strncpy(mirror_types[2],types1[0],5);
856 strncpy(mirror_types[3],types1[2],5);
857 k = find_match(4,mirror_types,item,&backwards);
858 if (k >= 0) return k;
859
860 /* d b c a */
861
862 *rearrange_ptr = 3;
863 strncpy(mirror_types[2],types1[2],5);
864 strncpy(mirror_types[3],types1[0],5);
865 k = find_match(4,mirror_types,item,&backwards);
866 if (k >= 0) return k;
867
868 /* c b a d */
869
870 *rearrange_ptr = 4;
871 strncpy(mirror_types[0],types1[2],5);
872 strncpy(mirror_types[2],types1[0],5);
873 strncpy(mirror_types[3],types1[3],5);
874 k = find_match(4,mirror_types,item,&backwards);
875 if (k >= 0) return k;
876
877 /* c b d a */
878
879 *rearrange_ptr = 5;
880 strncpy(mirror_types[2],types1[3],5);
881 strncpy(mirror_types[3],types1[0],5);
882 k = find_match(4,mirror_types,item,&backwards);
883 return k;
884 }
885
rearrange_improper(int ooptype,int rearrange)886 void rearrange_improper(int ooptype,int rearrange)
887 {
888 int i,j,temp[4];
889
890 for (i=0; i < 4; i++) temp[i] = ooptypes[ooptype].types[i];
891
892 switch (rearrange) {
893 case 1:
894 ooptypes[ooptype].types[0] = temp[0];
895 ooptypes[ooptype].types[2] = temp[3];
896 ooptypes[ooptype].types[3] = temp[2];
897 for (i=0; i < total_no_oops; i++) {
898 if (oops[i].type == ooptype) {
899 for (j=0; j < 4; j++) temp[j] = oops[i].members[j];
900 oops[i].members[2] = temp[3];
901 oops[i].members[3] = temp[2];
902 }
903 }
904 break;
905 case 2:
906 ooptypes[ooptype].types[0] = temp[3];
907 ooptypes[ooptype].types[2] = temp[0];
908 ooptypes[ooptype].types[3] = temp[2];
909 for (i=0; i < total_no_oops; i++) {
910 if (oops[i].type == ooptype) {
911 for (j=0; j < 4; j++) temp[j] = oops[i].members[j];
912 oops[i].members[0] = temp[3];
913 oops[i].members[2] = temp[0];
914 oops[i].members[3] = temp[2];
915 }
916 }
917 break;
918 case 3:
919 ooptypes[ooptype].types[0] = temp[3];
920 ooptypes[ooptype].types[2] = temp[2];
921 ooptypes[ooptype].types[3] = temp[0];
922 for (i=0; i < total_no_oops; i++) {
923 if (oops[i].type == ooptype) {
924 for (j=0; j < 4; j++) temp[j] = oops[i].members[j];
925 oops[i].members[0] = temp[3];
926 oops[i].members[2] = temp[2];
927 oops[i].members[3] = temp[0];
928 }
929 }
930 break;
931 case 4:
932 ooptypes[ooptype].types[0] = temp[2];
933 ooptypes[ooptype].types[2] = temp[0];
934 ooptypes[ooptype].types[3] = temp[3];
935 for (i=0; i < total_no_oops; i++) {
936 if (oops[i].type == ooptype) {
937 for (j=0; j < 4; j++) temp[j] = oops[i].members[j];
938 oops[i].members[0] = temp[2];
939 oops[i].members[2] = temp[0];
940 oops[i].members[3] = temp[3];
941 }
942 }
943 break;
944 case 5:
945 ooptypes[ooptype].types[0] = temp[2];
946 ooptypes[ooptype].types[2] = temp[3];
947 ooptypes[ooptype].types[3] = temp[0];
948 for (i=0; i < total_no_oops; i++) {
949 if (oops[i].type == ooptype) {
950 for (j=0; j < 4; j++) temp[j] = oops[i].members[j];
951 oops[i].members[0] = temp[2];
952 oops[i].members[2] = temp[3];
953 oops[i].members[3] = temp[0];
954 }
955 }
956 break;
957 default:
958 break;
959 }
960 }
961
find_trigonal_body_data(char types1[][5],struct FrcFieldItem item)962 int find_trigonal_body_data(char types1[][5],struct FrcFieldItem item)
963 {
964 int k,backwards;
965 char mirror_types[4][5];
966
967 backwards = -1;
968
969 /* a b c d */
970
971 k = find_match(4,types1,item,&backwards);
972 if (k >= 0) return k;
973
974 /* a b d c */
975
976 strncpy(mirror_types[0],types1[0],5);
977 strncpy(mirror_types[1],types1[1],5);
978 strncpy(mirror_types[2],types1[3],5);
979 strncpy(mirror_types[3],types1[2],5);
980 k = find_match(4,mirror_types,item,&backwards);
981 if (k >= 0) return k;
982
983 /* d b a c */
984
985 strncpy(mirror_types[0],types1[3],5);
986 strncpy(mirror_types[2],types1[0],5);
987 strncpy(mirror_types[3],types1[2],5);
988 k = find_match(4,mirror_types,item,&backwards);
989 if (k >= 0) return k;
990
991 /* d b c a */
992
993 strncpy(mirror_types[2],types1[2],5);
994 strncpy(mirror_types[3],types1[0],5);
995 k = find_match(4,mirror_types,item,&backwards);
996 if (k >= 0) return k;
997 /* c b a d */
998
999 strncpy(mirror_types[0],types1[2],5);
1000 strncpy(mirror_types[2],types1[0],5);
1001 strncpy(mirror_types[3],types1[3],5);
1002 k = find_match(4,mirror_types,item,&backwards);
1003 if (k >= 0) return k;
1004
1005 /* c b d a */
1006
1007 strncpy(mirror_types[2],types1[3],5);
1008 strncpy(mirror_types[3],types1[0],5);
1009 k = find_match(4,mirror_types,item,&backwards);
1010 return k;
1011 }
1012
find_angleangle_data(char types1[][5],struct FrcFieldItem item,int kloc[3])1013 int find_angleangle_data(char types1[][5],struct FrcFieldItem item,int kloc[3])
1014 {
1015 int k,backwards = -1;
1016 char mirror_types[4][5];
1017
1018 strncpy(mirror_types[1],types1[1],5);
1019
1020 /* go for first parameter a b c d or d b c a */
1021
1022 k = find_match(4,types1,item,&backwards);
1023 if (k < 0) {
1024 strncpy(mirror_types[0],types1[3],5);
1025 strncpy(mirror_types[2],types1[2],5);
1026 strncpy(mirror_types[3],types1[0],5);
1027 k = find_match(4,mirror_types,item,&backwards);
1028 }
1029 kloc[0] = k;
1030
1031 /* go for second parameter d b a c or c b a d */
1032
1033 strncpy(mirror_types[0],types1[3],5);
1034 strncpy(mirror_types[2],types1[0],5);
1035 strncpy(mirror_types[3],types1[2],5);
1036 k = find_match(4,mirror_types,item,&backwards);
1037 if (k < 0) {
1038 strncpy(mirror_types[0],types1[2],5);
1039 strncpy(mirror_types[3],types1[3],5);
1040 k = find_match(4,mirror_types,item,&backwards);
1041 }
1042 kloc[1] = k;
1043
1044 /* go for third parameter a b d c or c b d a */
1045
1046 strncpy(mirror_types[0],types1[0],5);
1047 strncpy(mirror_types[2],types1[3],5);
1048 strncpy(mirror_types[3],types1[2],5);
1049 k = find_match(4,mirror_types,item,&backwards);
1050 if (k < 0) {
1051 strncpy(mirror_types[0],types1[2],5);
1052 strncpy(mirror_types[3],types1[0],5);
1053 k = find_match(4,mirror_types,item,&backwards);
1054 }
1055 kloc[2] = k;
1056 k = 0;
1057 if ((kloc[0] < 0) && (kloc[1] < 0) && (kloc[2] < 0)) k = -1;
1058 return k;
1059 }
1060
find_match(int n,char types1[][5],struct FrcFieldItem item,int * backwards_ptr)1061 int find_match(int n, char types1[][5],struct FrcFieldItem item,int
1062 *backwards_ptr)
1063 {
1064 int k,match;
1065
1066 match = 0;
1067 k=0;
1068
1069 /* Try for an exact match (no wildcards) first */
1070
1071 while (!match && (k < item.entries)) {
1072 if (match_types(n, 0,types1,item.data[k].ff_types,backwards_ptr) == 1)
1073 match = 1;
1074 else
1075 k++;
1076 }
1077
1078 /* Try again - allow wildcard matching */
1079
1080 if (!match) {
1081 k=0;
1082 while (!match && (k < item.entries)) {
1083 if (match_types(n,1,types1,item.data[k].ff_types,backwards_ptr) == 1)
1084 match = 1;
1085 else
1086 k++;
1087 }
1088 }
1089 if (match) return k;
1090 else return -1;
1091 }
1092
match_types(int n,int wildcard,char types1[][5],char types2[][5],int * backwards_ptr)1093 int match_types(int n,int wildcard,char types1[][5],char types2[][5],
1094 int *backwards_ptr)
1095 {
1096 int k,match;
1097
1098 /* Routine to match short arrays of characters strings which contain
1099 atom potential types. The arrays range from 1 to 4 (VDW or equivalences,
1100 bond, angle, dihedrals or oops). There are potentially four ways the
1101 arrays can match: exact match (forwards), exact match when one array is
1102 run backwards (backwards), forwards with wildcard character match allowed
1103 (forwards *) and finally backwards with wildcard character match
1104 (backwards *). If the variable, backwards (pointed by backwards_ptr)
1105 is -1, then the backwards options are not to be used (such when
1106 matching oop types)
1107 */
1108
1109
1110 if (wildcard == 0) {
1111
1112 /* forwards */
1113
1114 k=0;
1115 match = 1;
1116 while (match && (k < n)) {
1117 if (strncmp(types1[k],types2[k],5) == 0)
1118 k++;
1119 else
1120 match = 0;
1121 }
1122 } else {
1123
1124 /* forwards * */
1125
1126 k=0;
1127
1128 match = 1;
1129 while (match && (k < n)) {
1130 if ((strncmp(types1[k],types2[k],5) == 0) ||
1131 (types2[k][0] == '*'))
1132 k++;
1133 else
1134 match = 0;
1135 }
1136 }
1137
1138 if (match) {
1139 *backwards_ptr = 0;
1140 return 1;
1141 }
1142 if ((n < 2) || (*backwards_ptr == -1)) return 0;
1143
1144 if (wildcard == 0) {
1145
1146 /* backwards */
1147
1148 k=0;
1149 match = 1;
1150 while (match && (k < n)) {
1151 if (strncmp(types1[n-k-1],types2[k],5) == 0)
1152 k++;
1153 else
1154 match = 0;
1155 }
1156 } else {
1157
1158 /* backwards * */
1159
1160 k=0;
1161 match = 1;
1162 while (match && (k < n)) {
1163 if ((strncmp(types1[n-k-1],types2[k],5) == 0) ||
1164 (types2[k][0] == '*') )
1165 k++;
1166 else
1167 match = 0;
1168 }
1169 }
1170
1171 if (match) {
1172 *backwards_ptr = 1;
1173 return 1;
1174 } else return 0;
1175 }
1176
get_r0(int typei,int typej)1177 double get_r0(int typei,int typej)
1178 {
1179 int k,match;
1180 double r;
1181
1182 k=0;
1183 match=0;
1184 r = 0.0;
1185
1186 while (!match && (k < no_bond_types)) {
1187 if (((typei == bondtypes[k].types[0]) &&
1188 (typej == bondtypes[k].types[1])) ||
1189 ((typej == bondtypes[k].types[0]) &&
1190 (typei == bondtypes[k].types[1])) ) {
1191 r = bondtypes[k].params[0];
1192 match = 1;
1193 } else k++;
1194 }
1195
1196 if (match == 0)
1197 printf(" Unable to find r0 for types %d %d\n",typei,typej);
1198 return r;
1199 }
1200
get_t0(int typei,int typej,int typek)1201 double get_t0(int typei,int typej,int typek)
1202 {
1203 int k,match;
1204 double theta;
1205
1206 k=0;
1207 match=0;
1208 theta = 0.0;
1209
1210 while (!match && (k < no_angle_types)) {
1211 if (((typei == angletypes[k].types[0]) &&
1212 (typej == angletypes[k].types[1]) &&
1213 (typek == angletypes[k].types[2])) ||
1214 ((typek == angletypes[k].types[0]) &&
1215 (typej == angletypes[k].types[1]) &&
1216 (typei == angletypes[k].types[2])) ) {
1217 theta = angletypes[k].params[0];
1218 match = 1;
1219 } else k++;
1220 }
1221
1222 if (match == 0)
1223 printf(" Unable to find t0 for types %d %d %d\n",
1224 typei,typej,typek);
1225 return theta;
1226 }
1227
quo_cp()1228 int quo_cp()
1229 {
1230 char cp[] = "cp ";
1231 int i,type,found;
1232
1233 i = 0;
1234 type = -1;
1235 found = 0;
1236
1237 while (!found && (i < no_atom_types)) {
1238 if (strncmp(atomtypes[i].potential,cp,2) == 0) {
1239 found = 1;
1240 type = i;
1241 } else i++;
1242 }
1243
1244 return type;
1245 }
1246
get_equivs(int ic,char potential_types[][5],char equiv_types[][5])1247 void get_equivs(int ic,char potential_types[][5],char equiv_types[][5])
1248 {
1249 int i,k;
1250 switch (ic) {
1251 case 1:
1252 k = find_equiv_type(potential_types[0]);
1253 if (k > -1) strncpy(equiv_types[0],equivalence.data[k].ff_types[1],5);
1254 break;
1255
1256 case 2:
1257 for (i=0; i < 2; i++) {
1258 k = find_equiv_type(potential_types[i]);
1259 if (k > -1) strncpy(equiv_types[i],equivalence.data[k].ff_types[2],5);
1260 }
1261 break;
1262 case 3:
1263 for (i=0; i < 3; i++) {
1264 k = find_equiv_type(potential_types[i]);
1265 if (k > -1) strncpy(equiv_types[i],equivalence.data[k].ff_types[3],5);
1266 }
1267 break;
1268 case 4:
1269 for (i=0; i < 4; i++) {
1270 k = find_equiv_type(potential_types[i]);
1271 if (k > -1) strncpy(equiv_types[i],equivalence.data[k].ff_types[4],5);
1272 }
1273 break;
1274
1275 case 5:
1276 for (i=0; i < 4; i++) {
1277 k = find_equiv_type(potential_types[i]);
1278 if (k > -1)
1279 strncpy(equiv_types[i],equivalence.data[k].ff_types[5],5);
1280 }
1281 break;
1282 default:
1283 printf(" Requesting equivalences of unsupported type: %d\n",ic);
1284 condexit(26);
1285 break;
1286 }
1287 return;
1288 }
1289
find_equiv_type(char potential_type[5])1290 int find_equiv_type(char potential_type[5])
1291 {
1292 int j,k,match;
1293
1294 j = -1;
1295 k = 0;
1296 match = 0;
1297
1298 while (!match && (k < equivalence.entries)) {
1299 if (strncmp(potential_type,
1300 equivalence.data[k].ff_types[0],5) == 0) {
1301 match = 1;
1302 j = k;
1303 } else {
1304 k++;
1305 }
1306 }
1307 if (j < 0)
1308 printf(" Unable to find equivalent type for %s\n",potential_type);
1309 return j;
1310 }
1311