NAMD
GromacsTopFile.C
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <math.h>
5 #include <algorithm>
6 #include "common.h"
7 #include "ResizeArray.h"
8 #include "InfoStream.h"
9 #include "GromacsTopFile.h"
10 #include "InfoStream.h"
11 
12 #define JOULES_PER_CALORIE 4.184
13 #define ANGSTROMS_PER_NM 10
14 
15 /* A GromacsTopFile represents the information stored in a GROMACS
16  topolgy file. This is an immutable type. */
17 
18 /* XXX warning: this code contains a few algorithms which run in a
19  time of O(n^3) or so in the size of the topology file, since I
20  haven't bothered to do any sorting. I don't think it matters much,
21  since it still manages to run on the biggest files in less than a
22  second, and nothing (memory or time) depends on the total number of
23  atoms in the simulation - all that matters is the number that have
24  to be individually defined. */
25 
26 /* GromacsTopFile initializes this to represent the data stored in the
27  file <filename>, or exits on error. */
28 #define LINESIZE 200
29 
30 /* modes */
31 #define UNKNOWN 0
32 #define ATOMS 1
33 #define MOLECULETYPE 2
34 #define MOLECULES 3
35 #define SYSTEM 4
36 #define BONDS 5
37 #define BONDTYPES 6
38 #define ANGLES 7
39 #define ATOMTYPES 8
40 #define ANGLETYPES 9
41 #define DIHEDRALS 10
42 #define DIHEDRALTYPES 11
43 #define DEFAULTS 12
44 #define NONBOND 13
45 #define PAIRS 14
46 #define EXCLUSIONS 15
47 #ifndef CODE_REDUNDANT
48 #define CODE_REDUNDANT 0
49 #endif
50 
51 #define MIN_DEBUG_LEVEL 3
52 #include "Debug.h"
53 
54 /* Initialize variables for exclusion calculation */
55 /* JLai */
58 int numExclusion = 0;
59 int numLJPair = 0;
60 int numGaussPair = 0;
62 
63 
64 
66  fudgeLJ = fudgeQQ = 1.0;
67  /* open the file */
68  FILE *f = fopen(filename,"r");
69  char buf[LINESIZE];
70  char modename[20];
71  int mode;
72  if(f==NULL) {
73  sprintf(buf,"Error opening file '%s'",filename);
74  NAMD_die(buf);
75  }
76 
77  /* really bad parser XXX probably just works on the files we
78  happen to have REWRITE THIS SOON. It should allow for \- line
79  continuations, check for errors in the file, etc. */
80  while(fgets(buf,LINESIZE-1,f)) {
81  char testchar;
82  int i,j;
83 
84  /* defaults */
85  int nbfunc, combrule;
86  char genpairs[20];
87 
88  /* atom buffers */
89  int num, resnum, chargegp, typenum;
90  char type[NAMESIZE+1], restype[NAMESIZE+1], atomname[NAMESIZE+1];
91  char particletype[NAMESIZE+1];
92  float charge, mass, c6, c12, junkf;
93 
94  /* moltype buffers */
95  int nrexcl;
96  char molname[LONGNAMESIZE+1];
97 
98  /* molInst buffers */
99  int copies;
100  MolInst *moleculeinstance;
101 
102  /* general atomset buffers */
103  int atomi, atomj, atomk, atoml;
104  char typea[NAMESIZE+1],typeb[NAMESIZE+1],
105  typec[NAMESIZE+1],typed[NAMESIZE+1];
106  const char *tmptypea,*tmptypeb,*tmptypec,*tmptyped;
107  int funct, index;
108  float c0,c1;
109 
110  /* bonds */
111  float b0,kB,th0,kth;
112 
113  /* dihedrals */
114  float c[6];
115  int mult=0;
116 
117  /* check for comments */
118  if(sscanf(buf," %c",&testchar)==1) {
119  if(testchar == ';') continue;
120  }
121  else { /* this is a blank line */
122  continue;
123  }
124 
125  /* check for a new mode */
126  if(sscanf(buf," [ %19[^] ] ]",modename)==1) {
127  /* switch the mode */
128  if(0==strcmp(modename,"atoms")) mode = ATOMS;
129  else if(0==strcmp(modename,"atomtypes")) mode = ATOMTYPES;
130  else if(0==strcmp(modename,"moleculetype")) mode = MOLECULETYPE;
131  else if(0==strcmp(modename,"molecules")) mode = MOLECULES;
132  else if(0==strcmp(modename,"system")) mode = SYSTEM;
133  else if(0==strcmp(modename,"bonds")) mode = BONDS;
134  else if(0==strcmp(modename,"bondtypes")) mode = BONDTYPES;
135  else if(0==strcmp(modename,"angles")) mode = ANGLES;
136  else if(0==strcmp(modename,"angletypes")) mode = ANGLETYPES;
137  else if(0==strcmp(modename,"dihedrals")) mode = DIHEDRALS;
138  else if(0==strcmp(modename,"dihedraltypes")) mode = DIHEDRALTYPES;
139  else if(0==strcmp(modename,"defaults")) mode = DEFAULTS;
140  else if(0==strcmp(modename,"nonbond_params")) mode = NONBOND;
141  // JLai
142  else if(0==strcmp(modename,"pairs")) mode = PAIRS;
143  else if(0==strcmp(modename,"exclusions")) mode = EXCLUSIONS;
144  else {
145  fprintf(stderr,"Warning: unknown mode %s\n",modename);
146  mode = UNKNOWN;
147  }
148 
149  continue;
150  }
151 
152  /* now do the appropriate thing based on the current mode */
153  switch(mode) {
154  case SYSTEM:
155  systemName = strdup(buf);
156  break;
157 
158  case DEFAULTS:
159  i = sscanf(buf," %d %d %20s %f %f",
160  &nbfunc,&combrule,genpairs,&fudgeLJ,&fudgeQQ);
161  if(i < 3) { // didn't get enough parameters
162  fprintf(stderr,"syntax error in DEFAULTS\n");
163  exit(1);
164  }
165  if(nbfunc != 1) { // I don't know how it works with nbfunc=2
166  fprintf(stderr,"Non-bonded function != 1 unsupported in DEFAULTS\n");
167  exit(1);
168  }
169  if(combrule != 1) { // same here
170  fprintf(stderr,"Combination rule != 1 unsupported in DEFAULTS\n");
171  exit(1);
172  }
173  if(0==strcmp(genpairs,"yes")) {
174  genPairs=1;
175  if(i!=5) {
176  fprintf(stderr,"syntax error in DEFAULTS\n");
177  exit(1);
178  }
179  // else fudgeLJ and fudgeQQ got written automatically
180  }
181  else genPairs=0;
182 
183  break;
184 
185  case NONBOND:
186  if(5 != sscanf(buf," %5s %5s %d %f %f",
187  typea, typeb, &funct, &c6, &c12)) {
188  fprintf(stderr,"Syntax error in NONBOND\n");
189  exit(1);
190  }
191  // convert kJ/mol*nm6 to kcal/mol*A6 and ..12 ..12
192  c6 = c6/JOULES_PER_CALORIE*1E6;
193  c12= c12/JOULES_PER_CALORIE*1E12;
194  vdwTable.addType(typea,typeb,c6,c12);
195  break;
196 
197  case BONDS:
198  i = sscanf(buf," %d %d %d %f %f",
199  &atomi,&atomj,&funct,&c0,&c1);
200  atomi--; // shift right away to a zero-indexing
201  atomj--;
202  if(i==3) {
203  tmptypea = genericMols[genericMols.size()-1]->getAtom(atomi)->getType();
204  tmptypeb = genericMols[genericMols.size()-1]->getAtom(atomj)->getType();
205  /* find the index and parameters */
206  index = bondTable.getParams(tmptypea, tmptypeb, funct, &b0, &kB);
207  if(index==-1) {
208  fprintf(stderr,"Required bondtype %s--%s (function %d) not found.\n",
209  tmptypea,tmptypeb,funct);
210  exit(1);
211  }
212  }
213  else if(i==5) {
214  /* first set the values of b0 and kB correctly */
215  b0 = c0*ANGSTROMS_PER_NM; /* convert nm to A */
216  if(funct==1) { /* harmonic potential */
217  /* convert kJ/nm2 to kcal/A2 and use E=kx2 instead of half that. */
218  kB = c1/JOULES_PER_CALORIE/100/2;
219  }
220  else if(funct==2) { /* special fourth-order potential */
221  /* convert to the normal harmonic constant and kJ/nm2 to kcal/A2 */
222  kB = 2*c1*c0*c0/JOULES_PER_CALORIE/100;
223  kB /= 2; /* use the NAMD system where E=kx2 */
224  }
225  else {
226  fprintf(stderr,"I don't know what funct=%d means in BONDS\n",funct);
227  exit(1);
228  }
229  /* look up the index */
230  index = bondTable.getIndex(b0,kB,funct);
231  }
232  else {
233  fprintf(stderr,"Syntax error in BONDS\n");
234  exit(1);
235  }
236 
237  genericMols[genericMols.size()-1]->addBond(atomi,atomj,index);
238  break;
239 
240  case BONDTYPES:
241  if(5 != sscanf(buf," %5s %5s %d %f %f",
242  typea,typeb,&funct,&c0,&c1)) {
243  fprintf(stderr,"Syntax error in BONDTYPES\n");
244  exit(1);
245  }
246 
247  /* first set the values of b0 and kB correctly */
248  b0 = c0*10; /* convert nm to A */
249  if(funct==1) { /* harmonic potential */
250  /* convert kJ/nm2 to kcal/A2 and use E=kx2 instead of half that. */
251  kB = c1/JOULES_PER_CALORIE/100/2;
252  }
253  else if(funct==2) { /* special fourth-order potential */
254  /* convert to the normal harmonic constant and kJ/nm2 to kcal/A2 */
255  kB = 2*c1*c0*c0/JOULES_PER_CALORIE/100;
256  kB /= 2; /* use the NAMD system where E=kx2 */
257  }
258  else {
259  fprintf(stderr,"I don't know what funct=%d means in BONDTYPES\n",funct);
260  exit(1);
261  }
262 
263  bondTable.addType(typea,typeb,b0,kB,funct);
264  break;
265 
266  case ANGLES:
267  i = sscanf(buf," %d %d %d %d %f %f",
268  &atomi,&atomj,&atomk,&funct,&c0,&c1);
269  atomi--; // shift right away to a zero-indexing
270  atomj--;
271  atomk--;
272  if(i == 4) { /* we have to look up the last two parameters */
273  tmptypea = genericMols[genericMols.size()-1]->getAtom(atomi)->getType();
274  tmptypeb = genericMols[genericMols.size()-1]->getAtom(atomj)->getType();
275  tmptypec = genericMols[genericMols.size()-1]->getAtom(atomk)->getType();
276  /* find the index and parameters */
277  index = angleTable.getParams(tmptypea, tmptypeb, tmptypec,
278  funct, &th0, &kth);
279  if(index==-1) {
280  fprintf(stderr,
281  "Required angletype %s--%s--%s (function %d) not found.\n",
282  tmptypea,tmptypeb,tmptypec,funct);
283  exit(1);
284  }
285  }
286  else if(i == 6) {
287  /* first set the values of th0 and kth correctly */
288  if(funct == 1) {
289  th0 = c0; /* both are in degrees */
290  kth = c1/JOULES_PER_CALORIE/2; /* convert kJ/rad2 to kcal/rad2 and use E=kx2 */
291  }
292  else if(funct == 2) {
293  th0 = c0; /* both are in degrees */
294  /* convert G96 kJ to kcal/rad2 and use E=kx2 */
295  kth = sin(th0*PI/180)*sin(th0*PI/180)*c1/JOULES_PER_CALORIE/2;
296  }
297  else {
298  fprintf(stderr,"I don't know what funct=%d means in ANGLES\n",funct);
299  exit(1);
300  }
301  /* add the angle type to our table */
302  index = angleTable.getIndex(th0,kth,funct);
303  }
304  else {
305  fprintf(stderr,"Syntax error (%d args) in ANGLES: %s\n",i,buf);
306  exit(1);
307  }
308 
309  /* add the angle to our table */
310  genericMols[genericMols.size()-1]->addAngle(atomi,atomj,atomk,index);
311  break;
312 
313  case ANGLETYPES:
314  if(6 != sscanf(buf," %5s %5s %5s %d %f %f",
315  typea,typeb,typec,&funct,&c0,&c1)) {
316  fprintf(stderr,"Syntax error in ANGLETYPES\n");
317  exit(1);
318  }
319  /* first set the values of th0 and kth correctly */
320  if(funct == 1) {
321  th0 = c0; /* both are in degrees */
322  kth = c1/JOULES_PER_CALORIE/2; /* convert kJ/rad2 to kcal/rad2 and use E=kx2 */
323  }
324  else if(funct == 2) {
325  th0 = c0; /* both are in degrees */
326  /* convert G96 kJ to kcal/rad2 and use E=kx2 */
327  kth = sin(th0*PI/180)*sin(th0*PI/180)*c1/JOULES_PER_CALORIE/2;
328  }
329  else {
330  fprintf(stderr,"I don't know what funct=%d means in ANGLETYPES\n",
331  funct);
332  exit(1);
333  }
334  angleTable.addType(typea,typeb,typec,th0,kth,funct);
335  break;
336 
337  case DIHEDRALS:
338  i = sscanf(buf," %d %d %d %d %d %f %f %f %f %f %f",
339  &atomi,&atomj,&atomk,&atoml,&funct,
340  &c[0],&c[1],&c[2],&c[3],&c[4],&c[5]);
341  atomi--; // shift right away to a zero-indexing
342  atomj--;
343  atomk--;
344  atoml--;
345  if(i==5) { /* we have to look up the parameters */
346  tmptypea = genericMols[genericMols.size()-1]->getAtom(atomi)->getType();
347  tmptypeb = genericMols[genericMols.size()-1]->getAtom(atomj)->getType();
348  tmptypec = genericMols[genericMols.size()-1]->getAtom(atomk)->getType();
349  tmptyped = genericMols[genericMols.size()-1]->getAtom(atoml)->getType();
350  /* find the index and parameters */
351  index = dihedralTable.getParams(tmptypea, tmptypeb, tmptypec,
352  tmptyped, funct, c, &mult);
353  if(index==-1) {
354  fprintf(stderr,
355  "Required dihedraltype %s--%s--%s--%s (function %d) not found.\n",
356  tmptypea,tmptypeb,tmptypec,tmptyped,funct);
357  exit(1);
358  }
359  }
360  else if(i==7 || i==8 || i==11) { /* the parameters are given */
361  if(funct==1 || funct==2) { /* we should have two parameters */
362  if(i!=7+(funct==1)) { /* plus a multiplicity for funct==1 */
363  fprintf(stderr,"Must have 7 args for funct=1,2\n");
364  exit(1);
365  }
366  c[0] = c[0]; /* both in deg */
367  if(i==7) {
368  c[1] = c[1]/(2*JOULES_PER_CALORIE); /* convert kJ to kcal and still use E=kx2*/
369  } else if (i==8 || i==11) {
370  c[1] = c[1]/(1*JOULES_PER_CALORIE); /* convert kJ to kcal and still use E=kx2*/
371  }
372  //c[1] = c[1]/(2*JOULES_PER_CALORIE); /* convert kJ to kcal and still use E=kx2*/
373  /* for funct==1 these are both divided by rad^2 */
374  if(funct==1) {
375  mult=(int)(c[2]+0.5); /* round to nearest integer :p */
376  }
377  }
378  else if(funct==3) {
379  if(i!=11){
380  fprintf(stderr,"Must have 11 args for funct=3\n");
381  exit(1);
382  }
383 
384  for(j=0;j<6;j++) {
385  c[j] = c[j]/JOULES_PER_CALORIE; /* convert kJ to kcal and E=Cn cos^n */
386  }
387  }
388  else {
389  fprintf(stderr,
390  "I don't know what funct=%d means in DIHEDRALS\n",funct);
391  exit(1);
392  }
393  index = dihedralTable.getIndex(c,mult,funct);
394  }
395  else {
396  fprintf(stderr,"Syntax error (%d args) in DIHEDRALS: %s\n",i,buf);
397  exit(1);
398  }
399 
400  /* add the dihedrals to our table */
401  genericMols[genericMols.size()-1]->addDihedral(atomi,atomj,atomk,atoml,
402  index);
403  break;
404  case DIHEDRALTYPES:
405  i = sscanf(buf," %5s %5s %d %f %f %f %f %f %f",
406  typea,typeb,&funct,&c[0],&c[1],&c[2],&c[3],&c[4],&c[5]);
407  if(funct == 1 || funct == 2) {
408  if(i!=5+(funct==1)) { /* 6 for f=2, 5 for f=1 */
409  fprintf(stderr,"Syntax error in DIHEDRALTYPES: %s\n",buf);
410  exit(1);
411  }
412  c[0] = c[0]; /* both in deg */
413  c[1] = c[1]/JOULES_PER_CALORIE; /* convert kJ to kcal and still use E=kx2*/
414  /* for funct==1 these are both divided by rad^2 */
415  if(funct==1) {
416  mult=(int)(c[2]+0.5); /* round to nearest integer :p */
417  }
418  }
419  else if(funct == 3) {
420  if(i!=9) {
421  fprintf(stderr,"Syntax error in DIHEDRALTYPES\n");
422  exit(1);
423  }
424  for(j=0;j<6;j++) {
425  c[j] = c[j]/JOULES_PER_CALORIE; /* convert kJ to kcal and E=Cn cos^n */
426  }
427  }
428  else {
429  fprintf(stderr,"I don't know what funct=%d means in DIHEDRALTYPES\n",
430  funct);
431  exit(1);
432  }
433  dihedralTable.addType(typea,typeb,c,mult,funct);
434  break;
435 
436  case ATOMS:
437  i = sscanf(buf," %d %5s %d %5s %5s %d %f %f",
438  &num, type, &resnum, restype,
439  atomname, &chargegp, &charge, &mass);
440  if(i==7) { /* XXX temporary - I should be able to get more
441  params */
442  typenum = atomTable.getParams(type,&mass,&junkf,&junkf,&junkf);
443  i=8;
444  }
445  else {
446  if(i!=8) {
447  fprintf(stderr,"Syntax error in ATOMS\n");
448  exit(1);
449  }
450  // just get the type number
451  typenum = atomTable.getParams(type,&junkf,&junkf,&junkf,&junkf);
452  }
453  genericMols[genericMols.size()-1]->addAtom(type,typenum,resnum,
454  restype,atomname,charge,mass);
455  break;
456 
457  case ATOMTYPES:
458  if(6 != sscanf(buf," %5s %f %f %5s %f %f",type,&mass,&charge,
459  particletype,&c6,&c12)) {
460  fprintf(stderr,"Syntax error in ATOMTYPES: %s\n",buf);
461  exit(1);
462  }
463  /* conversions:
464  c6 - kJ/mol nm6 -> kcal/mol A6
465  c12 - kJ/mol nm12 -> kcal/mol A12 */
466  atomTable.addType(type,mass,charge,
467  c6/(JOULES_PER_CALORIE)*1E6,
468  c12/(JOULES_PER_CALORIE)*1E12);
469  break;
470 
471  case MOLECULETYPE:
472  if(2!=sscanf(buf," %20s %d",molname,&nrexcl)) {
473  fprintf(stderr,"Syntax error in MOLECULETYPE: %s\n",buf);
474  exit(1);
475  }
476 
477  /* add another generic molecule holder */
478  genericMols.add(new GenericMol(molname));
479  break;
480 
481  case MOLECULES:
482  if(2!=sscanf(buf," %20s %d",molname,&copies)) {
483  fprintf(stderr,"Syntax error in MOLECULES: %s\n",buf);
484  exit(1);
485  }
486 
487  /* search for the specified molecule and make a molInst of it*/
488  moleculeinstance = NULL;
489  for(i=0;i<genericMols.size();i++) {
490  if(0==strcmp(molname,genericMols[i]->getName())) {
491  moleculeinstance = new MolInst(genericMols[i],copies);
492  break;
493  }
494  }
495 
496  if(moleculeinstance==NULL) {
497  fprintf(stderr,"Molecule %s undefined.\n",molname);
498  exit(1);
499  }
500 
501  /* put it at the end of the list */
502  molInsts.add(moleculeinstance);
503 
504  break;
505  case PAIRS:
506  int indexA;
507  int indexB;
508  int pairFunction;
509  Real fA;
510  Real fB;
511  Real fC;
512  Real fD;
513  Real fE;
514  Real fF;
515 
516  int countVariables;
517  countVariables = sscanf(buf," %d %d %d %f %f %f %f %f %f",&indexA,&indexB,&pairFunction,&fA,&fB,&fC,&fD,&fE,&fF);
518 
519  if ((countVariables >= 4 && countVariables >= 10)) {
520  fprintf(stderr,"Syntax error in PAIRS: %s\n",buf);
521  exit(1);
522  }
523 
524  // Shift the atom indices to be zero-based
525  indexA--;
526  indexB--;
527 
528  // Make sure that the indexA is less than indexB
529  /*if ( indexA > indexB ) {
530  int tmpIndex = indexA;
531  indexB = indexA;
532  indexA = tmpIndex;
533  }*/
534 
535  if (pairFunction == 1) {
536 
537  // LJ code
538  fA = (fA/JOULES_PER_CALORIE)*1E6;
539  fB= (fB/JOULES_PER_CALORIE)*1E12;
540  pairTable.addPairLJType2(indexA,indexB,fA,fB);
541  } else if (pairFunction == 5){
542 
543  // Bare Gaussian potential
544  fA = (fA/JOULES_PER_CALORIE); //-->gaussA
545  fB = (fB*ANGSTROMS_PER_NM); //-->gaussMu1
546  if(fC == 0) {
547  char buff[100];
548  sprintf(buff,"GromacsTopFile.C::Attempting to load zero into gaussSigma. Please check the pair: %s\n",buf);
549  NAMD_die(buff);
550  }
551  if(fC < 0 && !bool_negative_number_warning_flag) {
552  iout << iWARN << "Attempting to load a negative standard deviation into the gaussSigma. Taking the absolute value of the standard deviation.";
554  }
555  fC = (fC*ANGSTROMS_PER_NM); //-->gaussSigma1
556  fC = 1.0/(2 * fC * fC); // Normalizes sigma
557  pairTable.addPairGaussType2(indexA,indexB,fA,fB,fC);
558  } else if (pairFunction == 6) {
559 
560  // Combined Gaussian + repulsive r^-12 potential
561  fA = (fA/JOULES_PER_CALORIE); //-->gaussA
562  fB = (fB*ANGSTROMS_PER_NM); //-->gaussMu1
563  if(fC == 0) {
564  char buff[100];
565  sprintf(buff,"GromacsTopFile.C::Attempting to load zero into gaussSigma. Please check the pair: %s\n",buf);
566  NAMD_die(buff);
567  }
568  if(fC < 0 && !bool_negative_number_warning_flag) {
569  iout << iWARN << "Attempting to load a negative standard deviation into the gaussSigma. Taking the absolute value of the standard deviation.";
571  }
572  fC = (fC*ANGSTROMS_PER_NM); //-->gaussSigma1
573  fC = 1.0/(2 * fC * fC); // Normalizes sigma
574  fD = (fD*ANGSTROMS_PER_NM); //-->gaussRepulsive
575  pairTable.addPairGaussType2(indexA,indexB,fA,fB,fC,fD);
576  } else if (pairFunction == 7) {
577 
578  // Double well Guassian function
579  fA = (fA/JOULES_PER_CALORIE); //-->gaussA
580  fB = (fB*ANGSTROMS_PER_NM); //-->gaussMu1
581  if(fC == 0 || fE == 0) {
582  char buff[100];
583  sprintf(buff,"GromacsTopFile.C::Attempting to load zero into gaussSigma. Please check the pair: %s\n",buf);
584  NAMD_die(buff);
585  }
586  if((fC < 0 || fE < 0)&& !bool_negative_number_warning_flag) {
587  iout << iWARN << "Attempting to load a negative standard deviation into the gaussSigma. Taking the absolute value of the standard deviation.";
589  }
590  fC = (fC*ANGSTROMS_PER_NM); //-->gaussSigma1
591  fC = 1.0/(2 * fC * fC); // Normalizes sigma
592  fD = (fD*ANGSTROMS_PER_NM); //-->gaussMu2
593  fE = (fE*ANGSTROMS_PER_NM); //-->gaussSigma2
594  fE = 1.0/(2 * fE * fE); // Normalizes sigma
595  fF = (fE*ANGSTROMS_PER_NM); //-->gaussRepulsive
596  pairTable.addPairGaussType2(indexA,indexB,fA,fB,fC,fD,fE,fF);
597  } else {
598 
599  // Generic error statement
600  fprintf(stderr,"Unknown pairFunction in GromacsTopFile.C under the PAIRS section: %d\n",pairFunction);
601  }
602  break;
603  case EXCLUSIONS:
604  /* Start of JLai modifications August 16th, 2012 */
605  if(2 != sscanf(buf," %d %d ",&atomi,&atomj)) {
606  fprintf(stderr,"Syntax error in EXCLUSIONS: %s\n",buf);
607  exit(1);
608  }
609 
610  // Shift the atom indices to be zero-based
611  atomi--;
612  atomj--;
613 
614  /*Load exclusion information into file*/
615  exclusions_atom_i.add(atomi);
616  exclusions_atom_j.add(atomj);
617  numExclusion++;
618 
619  /* Reading in exclusions information from file and loading */
620  break;
621  // End of JLai modifications August 16th, 2012
622  }
623  }
624 
625  fclose(f);
626 }
627 
628 /* returns the number of bonds in the file */
630  int n=0,i;
631  for(i=0;i<molInsts.size();i++) {
632  n += molInsts[i]->getNum() *
633  molInsts[i]->getMol()->getNumBonds();
634  }
635  return n;
636 }
637 
638 /* returns the number of angles in the file */
640  int n=0,i;
641  for(i=0;i<molInsts.size();i++) {
642  n += molInsts[i]->getNum() *
643  molInsts[i]->getMol()->getNumAngles();
644  }
645  return n;
646 }
647 
648 /* returns the number of dihedral angles in the file */
650  int n=0,i;
651  for(i=0;i<molInsts.size();i++) {
652  n += molInsts[i]->getNum() *
653  molInsts[i]->getMol()->getNumDihedrals();
654  }
655  return n;
656 }
657 
658 /* JLai -- August 16th, 2012 modifications*/
659 /* returns the number of pair interactions in the file */
661  int numPair = 0;
662  numPair = numLJPair + numGaussPair;
663  return numPair;
664 }
665 
667  return numLJPair;
668 }
669 
671  return numGaussPair;
672 }
673 
674 /* return the number of exclusions in the file */
676  return numExclusion;
677 }
678 
679 /* return the list of exclusions from the file */
680 void GromacsTopFile::getExclusions(int* atomi, int* atomj) const {
681  for(int i =0; i < exclusions_atom_i.size(); i++) {
682  atomi[i] = exclusions_atom_i[i];
683  atomj[i] = exclusions_atom_j[i];
684  }
685  return;
686 }
687 /* End of JLai modifications */
688 
689 /* getBond puts the information about bond number <num> into the
690  spaces pointed to by the other arguments. Bond number 0 is the
691  first bond in the list.
692  For atomi and atomj, 1 is the first atom in the list. */
693 void GromacsTopFile::getBond(int n, int *atomi, int *atomj, int *bondtype) const {
694  /* figure out what instance this bond is in */
695  int nbonds=0;
696  int natoms=0;
697  int i;
698  for(i=0;i<molInsts.size();i++) {
699  int molbonds = molInsts[i]->getMol()->getNumBonds();
700  int molatoms = molInsts[i]->getMol()->getNumAtoms();
701  int newbonds = molInsts[i]->getNumBonds();
702  int newatoms = molInsts[i]->getNumAtoms();
703 
704  if(nbonds+newbonds-1 >= n) {
705  /* the bond is in this MolInst */
706  int localnum = (n-nbonds) % molbonds; /* number within this inst */
707  int instnum = (n-nbonds) / molbonds; /* number of instances before */
708  int addatoms = natoms+instnum*molatoms; /* extra atoms to add to atomi */
709 
710  const GenericBond *b = molInsts[i]->getMol()->getBond(localnum);
711 
712  *atomi = b->getAtomi() + addatoms;
713  *atomj = b->getAtomj() + addatoms;
714  *bondtype = b->getType();
715  break;
716  }
717 
718  nbonds += newbonds;
719  natoms += newatoms;
720  }
721 }
722 
723 /* getAngle puts the information about angle number <num> into the
724  spaces pointed to by the other arguments. Angle number 0 is the
725  first angle in the list.
726 */
727 void GromacsTopFile::getAngle(int n, int *atomi, int *atomj, int *atomk,
728  int *angletype) const {
729  /* figure out what instance this angle is in */
730  int nangles=0;
731  int natoms=0;
732  int i;
733  for(i=0;i<molInsts.size();i++) {
734  int molangles = molInsts[i]->getMol()->getNumAngles();
735  int molatoms = molInsts[i]->getMol()->getNumAtoms();
736  int newangles = molInsts[i]->getNumAngles();
737  int newatoms = molInsts[i]->getNumAtoms();
738 
739  if(nangles+newangles-1 >= n) {
740  /* the angle is in this MolInst */
741  int localnum = (n-nangles) % molangles; /* number within this inst */
742  int instnum = (n-nangles) / molangles; /* number of instances before */
743  int addatoms = natoms+instnum*molatoms; /* extra atms to add to atmi */
744 
745  const GenericAngle *a = molInsts[i]->getMol()->getAngle(localnum);
746 
747  *atomi = a->getAtomi() + addatoms;
748  *atomj = a->getAtomj() + addatoms;
749  *atomk = a->getAtomk() + addatoms;
750  *angletype = a->getType();
751  break;
752  }
753 
754  nangles += newangles;
755  natoms += newatoms;
756  }
757 }
758 
759 /* getDihedral puts the information about dihedral number <num> into
760  the spaces pointed to by the other arguments. Dihedral number 0
761  is the first angle in the list. */
762 void GromacsTopFile::getDihedral(int n, int *atomi, int *atomj, int *atomk,
763  int *atoml, int *type) const {
764  /* figure out what instance this angle is in */
765  int ndihedrals=0;
766  int natoms=0;
767  int i;
768  for(i=0;i<molInsts.size();i++) {
769  int moldihedrals = molInsts[i]->getMol()->getNumDihedrals();
770  int molatoms = molInsts[i]->getMol()->getNumAtoms();
771  int newdihedrals = molInsts[i]->getNumDihedrals();
772  int newatoms = molInsts[i]->getNumAtoms();
773 
774  if(ndihedrals+newdihedrals-1 >= n) {
775  /* the dihedral is in this MolInst */
776  int localnum = (n-ndihedrals) % moldihedrals; /* number in this inst */
777  int instnum = (n-ndihedrals) / moldihedrals; /* number of insts before */
778  int addatoms = natoms+instnum*molatoms; /*extra atms to add to atmi*/
779 
780  const GenericDihedral *a = molInsts[i]->getMol()->getDihedral(localnum);
781 
782  *atomi = a->getAtomi() + addatoms;
783  *atomj = a->getAtomj() + addatoms;
784  *atomk = a->getAtomk() + addatoms;
785  *atoml = a->getAtoml() + addatoms;
786  *type = a->getType();
787  break;
788  }
789 
790  ndihedrals += newdihedrals;
791  natoms += newatoms;
792  }
793 }
794 
795 /* getNumAtoms returns the number of atoms stored in the file. */
797  int n=0;
798  int i;
799  for(i=0;i<molInsts.size();i++) {
800  n += molInsts[i]->getNum() *
801  molInsts[i]->getMol()->getNumAtoms();
802  }
803  return n;
804 }
805 
806 /* getAtom puts the information about atom number <n> into the
807  spaces pointed to by the other arguments. The string buffers
808  must be at least 11 characters long. */
810  int *residue_number, char *residue_name,
811  char *atom_name, char *atom_type, int *atom_typenum,
812  Real *charge, Real *mass)
813  const {
814  int natoms=0,n=num; // zero-indexed array
815  int resnum=0;
816  /* figure out what instance the atom is in, and what residue number
817  it has */
818  int i;
819 
820  for(i=0;i<molInsts.size();i++) {
821  int numnew = molInsts[i]->getNumAtoms(); /* atoms in this MolInst */
822  int resmol = molInsts[i]->getMol()->getNumRes(); /* # res/mol */
823  int newres = molInsts[i]->getNumRes(); /* residues in this MolInst */
824 
825  if(natoms+numnew-1 >= n) {
826  /* the atom is in this molInst */
827  int localnum = (n-natoms) % molInsts[i]->getMol()->getNumAtoms();
828  int instnum = (n-natoms) / molInsts[i]->getMol()->getNumAtoms();
829 
830  // getAtom is zero-indexed
831  const GenericAtom *a = molInsts[i]->getMol()->getAtom(localnum);
832  int residue = resnum + resmol*instnum + a->getResNum();
833 
834  *residue_number = residue;
835  strncpy(residue_name,a->getResType(),11);
836  strncpy(atom_name,a->getAtomName(),11);
837  strncpy(atom_type,a->getType(),11);
838  *charge=a->getCharge();
839  *mass=a->getMass();
840  *atom_typenum=a->getTypeNum();
841  break;
842  }
843 
844  /* go to the next molInst */
845  natoms += numnew;
846  resnum += newres;
847  }
848 }
849 
850 GenericAtom::GenericAtom(const char *theType, int theTypeNum, int theResNum,
851  const char *theResType,const char *theAtomName,
852  Real theCharge, Real theMass) {
853  strncpy(type,theType,NAMESIZE+1);
854  typeNum = theTypeNum;
855  resNum = theResNum;
856  strncpy(resType,theResType,NAMESIZE+1);
857  strncpy(atomName,theAtomName,NAMESIZE+1);
858  charge = theCharge;
859  mass = theMass;
860 }
861 
862 /* initializes this to be a bond between atoms <i> and <j> of type
863  <type>, with <next> pointing to the next bond in the list. */
864 GenericBond::GenericBond(int i, int j, int theType) {
865  atomi=i;
866  atomj=j;
867  type=theType;
868 }
869 
870 /* initializes this to be a angle between atoms <i>, <j>, and <k> of
871  type <type>, with <next> pointing to the next angle in the list. */
872 GenericAngle::GenericAngle(int i, int j, int k, int theType) {
873  atomi=i;
874  atomj=j;
875  atomk=k;
876  type=theType;
877 }
878 
879 /* initializes this to be a angle between atoms <i>, <j>, <k>, and
880  <l> of type <type> */
881 GenericDihedral::GenericDihedral(int i, int j, int k, int l, int theType) {
882  atomi=i;
883  atomj=j;
884  atomk=k;
885  atoml=l;
886  type=theType;
887 }
888 
889 /* adds a bond to the list */
890 void GenericMol::addBond(int atomi, int atomj, int type) {
891  bondList.add(new GenericBond(atomi, atomj, type));
892 }
893 
894 /* adds an angle to the list */
895 void GenericMol::addAngle(int atomi, int atomj, int atomk, int type) {
896  angleList.add(new GenericAngle(atomi, atomj, atomk, type));
897 }
898 
899 /* adds a dihedral to the list */
900 void GenericMol::addDihedral(int atomi, int atomj, int atomk,
901  int atoml, int type) {
902  dihedralList.add(new GenericDihedral(atomi, atomj, atomk,
903  atoml, type));
904 }
905 
906 /* adds an atom to the list */
907 void GenericMol::addAtom(const char *theType, int theTypeNum, int theResNum,
908  const char *theResType,const char *theAtomName,
909  Real theCharge, Real theMass) {
910 
911  atomList.add(new GenericAtom(theType,theTypeNum,theResNum,theResType,
912  theAtomName,theCharge,theMass));
913 }
914 
915 /* initializes this to be the molecule with name <name> */
916 GenericMol::GenericMol(const char *theName) {
917  name = strdup(theName);
918 }
919 
920 /* gets a bond */
921 const GenericBond *GenericMol::getBond(int n) const {
922  /* double-check */
923  if(n >= bondList.size() || n<0) {
924  fprintf(stderr,"Bond index %d out of bounds.\n",n);
925  exit(1);
926  }
927  return bondList[n];
928 }
929 
930 /* gets a angle */
931 const GenericAngle *GenericMol::getAngle(int n) const {
932  /* double-check */
933  if(n >= angleList.size() || n<0) {
934  fprintf(stderr,"Angle index %d out of bounds.\n",n);
935  exit(1);
936  }
937  return angleList[n];
938 }
939 
940 /* gets a dihedral */
942  /* double-check */
943  if(n >= dihedralList.size() || n<0) {
944  fprintf(stderr,"Dihedral index %d out of bounds.\n",n);
945  exit(1);
946  }
947  return dihedralList[n];
948 }
949 
950 /* gets an atom */
951 const GenericAtom *GenericMol::getAtom(int n) const {
952  /* double-check */
953  if(n >= atomList.size() || n<0) {
954  fprintf(stderr,"Atom index %d out of bounds for %s.\n",n,name);
955  exit(1);
956  }
957  return atomList[n];
958 }
959 
960 /* initializes this to represent <theNum> copies of <theMol> */
961 MolInst::MolInst(const GenericMol *theMol, int theNum) {
962  mol = theMol;
963  num = theNum;
964 }
965 
966 /* get the total number of various things here */
967 int MolInst::getNumAtoms() const {
968  return getMol()->getNumAtoms() * getNum();
969 }
970 int MolInst::getNumBonds() const {
971  return getMol()->getNumBonds() * getNum();
972 }
974  return getMol()->getNumAngles() * getNum();
975 }
977  return getMol()->getNumDihedrals() * getNum();
978 }
979 int MolInst::getNumRes() const {
980  return getMol()->getNumRes()
981  * getNum();
982 }
983 
984 /* this gets the index of a bond in the table (adding an entry if
985  none exists).
986  b0: the natural bond length in A.
987  kB: the spring constant in kcal/mol/A^2, where E=kx^2 to 1st order.
988  funct: 1 for harmonic potential, 2 for fourth-order GROMACS96
989  potential. */
990 int BondTable::getIndex(float b0, float kB, int funct) {
991  /* check to see if it is in the table already */
992  int i;
993  for(i=0;i<b0Array.size();i++) {
994  if(fabs(b0-b0Array[i])<0.00001 &&
995  fabs(kB-kBArray[i])<0.00001 &&
996  funct == functArray[i]) {
997  return i;
998  }
999  }
1000 
1001  /* nope, it wasn't in the table add a new element! */
1002  b0Array.add(b0);
1003  kBArray.add(kB);
1004  functArray.add(funct);
1005  typeaArray.add(NULL);
1006  typebArray.add(NULL);
1007  return b0Array.size()-1;
1008 }
1009 
1010 /* this gets the index of a angle in the table (adding an entry if
1011  none exists).
1012  b0: the natural angle length in A.
1013  kB: the spring constant in kcal/mol/A^2, where E=kx^2 to 1st order.
1014  funct: 1 for harmonic potential, 2 for fourth-order GROMACS96
1015  potential. */
1016 int AngleTable::getIndex(float th0, float kth, int funct) {
1017  /* check to see if it is in the table already */
1018  int i;
1019  for(i=0;i<th0Array.size();i++) {
1020  if(fabs(th0-th0Array[i])<0.00001 &&
1021  fabs(kth-kthArray[i])<0.00001 &&
1022  funct == functArray[i]) {
1023  return i;
1024  }
1025  }
1026 
1027  /* nope, it wasn't in the table add a new element! */
1028  th0Array.add(th0);
1029  kthArray.add(kth);
1030  functArray.add(funct);
1031  typeaArray.add(NULL);
1032  typebArray.add(NULL);
1033  typecArray.add(NULL);
1034  return th0Array.size()-1;
1035 }
1036 
1037 /* this function gets the index of a dihedral angle in the table
1038  (adding an entry if none exists). The required number of
1039  parameters (see notes on class DihedralTable) must be stored in
1040  the array <c> */
1041 int DihedralTable::getIndex(const float *c, int mult, int funct) {
1042  /* check to see if it is in the table already */
1043  int i,j,jmax;
1044  if(funct==1 || funct==2) { /* for these we only need two params */
1045  jmax=2;
1046  } else { /* for RB we need all 6 */
1047  jmax=6;
1048  }
1049 
1050  for(i=0;i<cArray.size();i++) {
1051  int mismatch=0;
1052  if(mult != multArray[i]) continue;
1053  if(funct != functArray[i]) continue;
1054 
1055  for(j=0;j<jmax;j++) {
1056  if(fabs(c[j]-cArray[i][j])>=0.00001) {
1057  mismatch=1;
1058  break;
1059  }
1060  }
1061  if(!mismatch) {
1062  /* all of the parameters matched */
1063  return i;
1064  }
1065  }
1066 
1067  /* nope, it wasn't in the table add a new element! */
1068  addType(NULL,NULL,c,mult,funct);
1069  return cArray.size()-1;
1070 }
1071 
1072 /* getBondParams puts the parameters for bond-type <num> into the
1073  spaces pointed to by the other arguments.
1074  b0 - natural length in A
1075  kB - spring constant in kcal/A^2 - E=kx^2
1076  funct - 1 for harmonic, 2 for special fourth-order GROMOS96 thing */
1077 void GromacsTopFile::getBondParams(int num, float *b0, float *kB, int *funct)
1078  const {
1079  bondTable.getParams(num,b0,kB,funct);
1080 }
1081 
1082 /* getAngleParams puts the parameters for angle-type <num> into the
1083  spaces pointed to by the other arguments. */
1084 void GromacsTopFile::getAngleParams(int num, float *th0, float *kth, int *funct)
1085  const {
1086  angleTable.getParams(num,th0,kth,funct);
1087 }
1088 
1089 void GromacsTopFile::getDihedralParams(int num, float *c, int *mult, int *funct)
1090  const {
1091  dihedralTable.getParams(num,c,mult,funct);
1092 }
1093 
1094 void GromacsTopFile::getPairLJArrays2(int *indexA, int *indexB, Real *pairC6, Real *pairC12) {
1095  pairTable.getPairLJArrays2(indexA, indexB, pairC6, pairC12);
1096 }
1097 
1098 void GromacsTopFile::getPairGaussArrays2(int *indexA, int *indexB, Real *gaussA, Real *gaussMu1,
1099  Real *gaussSigma1, Real *gaussMu2, Real *gaussSigma2,
1100  Real *gaussRepulsive) {
1101  pairTable.getPairGaussArrays2(indexA, indexB, gaussA, gaussMu1, gaussSigma1, gaussMu2, gaussSigma2, gaussRepulsive);
1102 }
1103 
1104 /* getParams puts the parameters for bond-type <num> into the
1105  spaces pointed to by the other arguments.
1106  b0 - natural length in A
1107  kB - spring constant in kcal/A^2 - E=kx^2
1108  funct - 1 for harmonic, 2 for special fourth-order GROMOS96 thing */
1109 void BondTable::getParams(int num, float *b0, float *kB, int *funct)
1110  const {
1111  *b0=b0Array[num];
1112  *kB=kBArray[num];
1113  *funct=functArray[num];
1114 }
1115 
1116 /* getParams puts the parameters for angle-type <num> into the
1117  spaces pointed to by the other arguments. */
1118 void AngleTable::getParams(int num, float *th0, float *kth, int *funct)
1119  const {
1120  *th0=th0Array[num];
1121  *kth=kthArray[num];
1122  *funct=functArray[num];
1123 }
1124 
1125 /* getParams puts the parameters for angle-type <num> into the
1126  spaces pointed to by the other arguments. The required number of
1127  parameters (see notes on class DihedralTable) will be stored in
1128  the array <c>, so make sure that c has size >= 6 */
1129 void DihedralTable::getParams(int num, float *c, int *mult, int *funct) const {
1130  int i;
1131 
1132  *funct=functArray[num]; /* first set the function */
1133 
1134  if(*funct==1 || *funct==2) { /* for these we only need two params */
1135  c[0]=cArray[num][0];
1136  c[1]=cArray[num][1];
1137  } else if(*funct==3) { /* for RB we need all 6 */
1138  for(i=0;i<6;i++) {
1139  c[i]=cArray[num][i];
1140  }
1141  } else { /* error */
1142  fprintf(stderr,"Bad function number %d - don't know what to do!\n",
1143  *funct);
1144  exit(1);
1145  }
1146 
1147  if(*funct==1) { /* return the multiplicity */
1148  *mult=multArray[num];
1149  }
1150 }
1151 
1152 /* this adds a entry for a angle type to the table, including the
1153  three atoms involved in the angle */
1154 void AngleTable::addType(const char *typea, const char *typeb,
1155  const char *typec, float th0, float kth, int funct) {
1156  typeaArray.add(strdup(typea));
1157  typebArray.add(strdup(typeb));
1158  typecArray.add(strdup(typec));
1159  th0Array.add(th0);
1160  kthArray.add(kth);
1161  functArray.add(funct);
1162 }
1163 
1164 /* this adds a entry for a angle type to the table, including the
1165  two atoms involved in the angle. The required number of
1166  parameters (see notes on class DihedralTable) must be stored in
1167  the array <c>. Note that these two angles really refer to either
1168  atoms A and D or B and C depending on the dihedral type. */
1169 void DihedralTable::addType(const char *typea, const char *typeb,
1170  const float *c, int mult, int funct) {
1171  float *cNew;
1172  int i;
1173 
1174  if(typea != NULL) typeaArray.add(strdup(typea));
1175  if(typeb != NULL) typebArray.add(strdup(typeb));
1176  functArray.add(funct);
1177 
1178  if(funct==1) multArray.add(mult);
1179  else multArray.add(0);
1180 
1181  if(funct==1 || funct==2) { /* for these we only need two params */
1182  cNew = new float[2];
1183  cNew[0]=c[0];
1184  cNew[1]=c[1];
1185  } else if(funct==3) { /* for RB we need all 6 */
1186  cNew = new float[6];
1187  for(i=0;i<6;i++) {
1188  cNew[i]=c[i];
1189  }
1190  } else { /* error */
1191  fprintf(stderr,"Bad function number %d - don't know what to do!\n",
1192  funct);
1193  exit(1);
1194  }
1195 
1196  /* now push the actual parameters */
1197  cArray.add(cNew);
1198 }
1199 
1200 /* This version looks up the angle by atom type - the direction of
1201  the types doesn't matter (it can be A--B--C or C--B--A)! If the
1202  specified angle/function combination cannot be found, it returns
1203  -1, otherwise it returns the index of the angletype. */
1204 int AngleTable::getParams(const char *typea, const char *typeb,
1205  const char *typec, int funct,
1206  float *th0, float *kth) const {
1207  int i;
1208  for(i=0;i<th0Array.size();i++) {
1209  if(typeaArray[i] == NULL || typebArray[i] == NULL
1210  || typecArray[i] == NULL) continue;
1211  if( (0==strcmp(typea,typeaArray[i]) && /* A--B--C */
1212  0==strcmp(typeb,typebArray[i]) &&
1213  0==strcmp(typec,typecArray[i])) /* or */
1214  || (0==strcmp(typec,typeaArray[i]) &&
1215  0==strcmp(typeb,typebArray[i]) && /* C--B--A */
1216  0==strcmp(typea,typecArray[i]))
1217  && funct==functArray[i]) {
1218  *th0 = th0Array[i];
1219  *kth = kthArray[i];
1220  return i;
1221  }
1222  }
1223  return -1;
1224 }
1225 
1226 /* see the (long) notes in the prototype definition */
1227 int DihedralTable::getParams(const char *typea, const char *typeb,
1228  const char *typec, const char *typed,
1229  int funct, float *c, int *mult) const {
1230  int i,j;
1231  const char *comparea, *compareb;
1232 
1233  if(funct == 1 || funct == 3) { /* for these, we use the inner atoms */
1234  comparea = typeb;
1235  compareb = typec;
1236  } else { /* use the outer atoms */
1237  comparea = typea;
1238  compareb = typed;
1239  }
1240 
1241  for(i=0;i<cArray.size();i++) {
1242  if(typeaArray[i] == NULL || typebArray[i] == NULL)
1243  continue; /* no atom types specified */
1244  if(functArray[i] != funct)
1245  continue; /* wrong function type */
1246 
1247  if( (0==strcmp(comparea,typeaArray[i]) && /* A--B */
1248  0==strcmp(compareb,typebArray[i])) /* or */
1249  || (0==strcmp(compareb,typeaArray[i]) &&
1250  0==strcmp(comparea,typebArray[i])) /* B--A */
1251  ) {
1252  if(funct==1 || funct==2) { /* for these we only need two params */
1253  c[0]=cArray[i][0];
1254  c[1]=cArray[i][1];
1255  if(funct==1) {
1256  *mult = multArray[i];
1257  }
1258  } else if(funct==3) { /* for RB we need all 6 */
1259  for(j=0;j<6;j++) {
1260  c[j]=cArray[i][j];
1261  }
1262  }
1263  return i;
1264  }
1265  }
1266  return -1;
1267 }
1268 
1269 /* this adds a entry for a bond type to the table, including the two
1270  atoms involved in the bond */
1271 void BondTable::addType(const char *typea, const char *typeb,
1272  float b0, float kB, int funct) {
1273  b0Array.add(b0);
1274  kBArray.add(kB);
1275  functArray.add(funct);
1276  typeaArray.add(strdup(typea));
1277  typebArray.add(strdup(typeb));
1278 }
1279 
1280 /* This version looks up the bond by atom type - the order of the
1281  types doesn't matter! */
1282 int BondTable::getParams(const char *typea, const char *typeb,
1283  int funct, float *b0, float *kB) const {
1284  int i;
1285  for(i=0;i<b0Array.size();i++) {
1286  if(typeaArray[i] == NULL || typebArray[i] == NULL) continue;
1287  if( (0==strcmp(typea,typeaArray[i]) &&
1288  0==strcmp(typeb,typebArray[i]))
1289  || (0==strcmp(typeb,typeaArray[i]) &&
1290  0==strcmp(typea,typebArray[i]))
1291  && funct==functArray[i]) {
1292  *b0 = b0Array[i];
1293  *kB = kBArray[i];
1294  return i;
1295  }
1296  }
1297  return -1;
1298 }
1299 
1300 /* this adds a entry for an atom type to the table */
1301 void AtomTable::addType(const char *type, float m, float q,
1302  float c6, float c12) {
1303  typeArray.add(strdup(type));
1304  mArray.add(m);
1305  qArray.add(q);
1306  c6Array.add(c6);
1307  c12Array.add(c12);
1308 }
1309 
1310 /* This looks up the atom type by number, returning it in the string
1311  <type>, which must be at least 11 characters long. */
1312 void AtomTable::getType(int num, char *type) const {
1313  if(num>=mArray.size() || num<0) {
1314  fprintf(stderr,"atomParam index %d out of bounds!\n",num);
1315  exit(1);
1316  }
1317  strncpy(type,typeArray[num],NAMESIZE+1);
1318 }
1319 
1320 /* This looks up the atom by type - if it cannot be found, this
1321  returns -1, otherwise this returns the index of the atomtype (for
1322  consistency - this is really a useless number.) */
1323 int AtomTable::getParams(const char *type, float *m, float *q,
1324  float *c6, float *c12) const {
1325  int i;
1326  for(i=0;i<mArray.size();i++) {
1327  if(typeArray[i]==NULL) {
1328  fprintf(stderr,"Found NULL atom type in list.\n");
1329  exit(1);
1330  }
1331  if(0==strcmp(typeArray[i],type)) {
1332  *m = mArray[i];
1333  *q = qArray[i];
1334  *c6 = c6Array[i];
1335  *c12 = c12Array[i];
1336  return i;
1337  }
1338  }
1339  return -1;
1340 }
1341 
1342 /* finds the index from the two interacting types, or returns -1 */
1343 int VDWTable::getIndex(const char *typea, const char *typeb) const {
1344  int i;
1345  for(i=0;i<c6Array.size();i++) {
1346  if((0==strcmp(typea,typeAArray[i]) &&
1347  0==strcmp(typeb,typeBArray[i])) ||
1348  (0==strcmp(typeb,typeAArray[i]) &&
1349  0==strcmp(typea,typeBArray[i]))) {
1350  return i;
1351  }
1352  }
1353  return -1;
1354 }
1355 
1356 /* these add VDW parameters to the list */
1357 void VDWTable::addType(const char *typea, const char *typeb, Real c6,
1358  Real c12) {
1359  int i;
1360 
1361  /* check to see if the pair is already in the table */
1362  i = getIndex(typea,typeb);
1363  if(i != -1) { /* it was in the table */
1364  c6Array[i] = c6;
1365  c12Array[i] = c12;
1366  }
1367  else { /* it wasn't in the table - add it! */
1368  typeAArray.add(strdup(typea));
1369  typeBArray.add(strdup(typeb));
1370  c6Array.add(c6);
1371  c12Array.add(c12);
1372  c6PairArray.add(0);
1373  c12PairArray.add(0);
1374  }
1375 }
1376 
1377 void VDWTable::add14Type(const char *typea, const char *typeb,
1378  Real c6pair, Real c12pair) {
1379  int i;
1380 
1381  /* check to see if the pair is already in the table */
1382  i = getIndex(typea,typeb);
1383  if(i != -1) { /* it was in the table */
1384  c6PairArray[i] = c6pair;
1385  c12PairArray[i] = c12pair;
1386  }
1387  else { /* it wasn't in the table - add it! */
1388  typeAArray.add(strdup(typea));
1389  typeBArray.add(strdup(typeb));
1390  c6PairArray.add(c6pair);
1391  c12PairArray.add(c12pair);
1392  c6Array.add(0);
1393  c12Array.add(0);
1394  }
1395 }
1396 
1397 /* this returns the VDW interaction parameters that were added to
1398  the list earlier, and returns the index or the parameters (a
1399  useless number) or -1 if it can't find the specified types. */
1400 int VDWTable::getParams(const char *typea, const char *typeb,
1401  Real *c6, Real *c12, Real *c6pair, Real *c12pair) const {
1402  int i;
1403  /* check to see if the pair is already in the table */
1404  i = getIndex(typea,typeb);
1405  if(i != -1) { /* it was in the table - return the parameters */
1406  *c6 = c6Array[i];
1407  *c12 = c12Array[i];
1408  *c6pair = c6PairArray[i];
1409  *c12pair = c12PairArray[i];
1410  }
1411  return i;
1412 }
1413 
1414 /* getVDWParams returs the Lennard-Jones bonding parameters for the
1415  specified two atom types, and the modified bonding parameters
1416  for 1-4 L-J interactons (c6pair, c12pair). */
1417 void GromacsTopFile::getVDWParams(int numa, int numb,
1418  Real *c6, Real *c12, Real *c6pair, Real *c12pair) const {
1419  int i,ret;
1420  Real c6a,c6b,c12a,c12b, mjunk, qjunk;
1421  char typea[11]="",typeb[11]="";
1422 
1423  /* get the string names corresponding to numa and numb */
1424  getAtomParams(numa,typea);
1425  getAtomParams(numb,typeb);
1426  if(typea[0]==0) { /* found a bug in my use of strncpy here once */
1427  fprintf(stderr,"Failed to get name of atom %d\n",numa);
1428  exit(1);
1429  }
1430  if(typeb[0]==0) {
1431  fprintf(stderr,"Failed to get name of atom %d\n",numb);
1432  exit(1);
1433  }
1434 
1435  /* first try - use the VDW table */
1436  i = vdwTable.getParams(typea,typeb,c6,c12,c6pair,c12pair);
1437 
1438  if(i==-1) {
1439  // QQ151069
1440  /*if ( !genPairs && numa != numb ) {
1441  iout << iWARN << "VDW table using combining rule for types "
1442  << typea << " " << typeb << "\n" << endi;
1443  }*/
1444 
1445  /* fallback - use the individual atom's parameters */
1446  ret = atomTable.getParams(typea, &mjunk, &qjunk, &c6a, &c12a);
1447  if(ret==-1) {
1448  fprintf(stderr,"Couldn't find atom type %s\n",typea);
1449  exit(1);
1450  }
1451  ret = atomTable.getParams(typeb, &mjunk, &qjunk, &c6b, &c12b);
1452  if(ret==-1) {
1453  fprintf(stderr,"Couldn't find atom type %s\n",typeb);
1454  exit(1);
1455  }
1456 
1457  *c6 = (float)(sqrt(c6a * c6b));
1458  *c12 = (float)(sqrt(c12a*c12b));
1459 
1460  /* this is the only reasonable option */
1461  *c6pair = *c6 * fudgeLJ;
1462  *c12pair = *c12 * fudgeLJ;
1463 
1464  }
1465 
1466 }
1467 
1468 /* Start of JLai modifications August 16th, 2012 */
1469 int PairTable::addPairLJType2(int indexA, int indexB, Real pairC6, Real pairC12) {
1470  GroLJPair glp;
1471  glp.indxLJA = indexA;
1472  glp.indxLJB = indexB;
1473  glp.c6pair = pairC6;
1474  glp.c12pair = pairC12;
1475  pairlistLJ.add(glp);
1476  numLJPair++;
1477  return 0;
1478 
1479  // Insert the second copy
1480  GroLJPair glp2;
1481  glp2.indxLJA = indexB;
1482  glp2.indxLJB = indexA;
1483  glp2.c6pair = pairC6;
1484  glp2.c12pair = pairC12;
1485  pairlistLJ.add(glp2);
1486  numLJPair++;
1487  return 0;
1488 }
1489 
1490 int PairTable::addPairGaussType2(int indexA, int indexB, Real gaussA, Real gaussMu1,
1491  Real gaussSigma1) {
1492  return addPairGaussType2(indexA, indexB, gaussA, gaussMu1, gaussSigma1, 0.0, 0.0, 0.0);
1493 }
1494 
1495 int PairTable::addPairGaussType2(int indexA, int indexB, Real gaussA, Real gaussMu1,
1496  Real gaussSigma1, Real gaussRepulsive) {
1497  return addPairGaussType2(indexA, indexB, gaussA, gaussMu1, gaussSigma1, 0.0, 0.0, gaussRepulsive);
1498 }
1499 
1500 int PairTable::addPairGaussType2(int indexA, int indexB, Real gaussA, Real gaussMu1,
1501  Real gaussSigma1, Real gaussMu2, Real gaussSigma2,
1502  Real gaussRepulsive) {
1503  GroGaussPair ggp;
1504  ggp.indxGaussA = indexA;
1505  ggp.indxGaussB = indexB;
1506  ggp.gA = gaussA;
1507  ggp.gMu1 = gaussMu1;
1508  ggp.giSigma1 = gaussSigma1;
1509  ggp.gMu2 = 0.0;
1510  ggp.giSigma2 = 0.0;
1511  ggp.gRepulsive = 0.0;
1512  pairlistGauss.add(ggp);
1513  numGaussPair++;
1514  GroGaussPair ggp2;
1515  ggp2.indxGaussA = indexB;
1516  ggp2.indxGaussB = indexA;
1517  ggp2.gA = gaussA;
1518  ggp2.gMu1 = gaussMu1;
1519  ggp2.giSigma1 = gaussSigma1;
1520  ggp2.gMu2 = 0.0;
1521  ggp2.giSigma2 = 0.0;
1522  ggp2.gRepulsive = 0.0;
1523  numGaussPair++;
1524  pairlistGauss.add(ggp2);
1525  return 0;
1526 }
1527 
1528 void PairTable::getPairLJArrays2(int *indexA, int *indexB, Real *pairC6, Real *pairC12) {
1529 
1530  std::sort(pairlistLJ.begin(),pairlistLJ.end(),GroLJCompare);
1532  for(int i = 0; i < numLJPair; i++) {
1533  indexA[i] = pairlistLJ[i].indxLJA;
1534  indexB[i] = pairlistLJ[i].indxLJB;
1535  pairC6[i] = pairlistLJ[i].c6pair;
1536  pairC12[i] = pairlistLJ[i].c12pair;
1537  }
1538 }
1539 
1540 void PairTable::getPairGaussArrays2(int *indexA, int *indexB, Real *gaussA, Real *gaussMu1,
1541  Real *gaussSigma1, Real *gaussMu2, Real *gaussSigma2,
1542  Real *gaussRepulsive) {
1543  std::sort(pairlistGauss.begin(),pairlistGauss.end(),GroGaussCompare);
1544  for(int i = 0; i < numGaussPair; i++) {
1545  indexA[i] = pairlistGauss[i].indxGaussA;
1546  indexB[i] = pairlistGauss[i].indxGaussB;
1547  gaussA[i] = pairlistGauss[i].gA;
1548  gaussMu1[i] = pairlistGauss[i].gMu1;
1549  gaussSigma1[i] = pairlistGauss[i].giSigma1;
1550  gaussMu2[i] = pairlistGauss[i].gMu2;
1551  gaussSigma2[i] = pairlistGauss[i].giSigma2;
1552  gaussRepulsive[i] = pairlistGauss[i].gRepulsive;
1553  }
1554 }
1555 
1557  if(A.indxLJA < B.indxLJA) {
1558  return true;
1559  } else if(A.indxLJA == B.indxLJA) {
1560  return (A.indxLJB < B.indxLJB);
1561  }
1562  return false;
1563 }
1564 
1566  if(A.indxGaussA < B.indxGaussA) {
1567  return true;
1568  } else if(A.indxGaussA == B.indxGaussA) {
1569  return (A.indxGaussB < B.indxGaussB);
1570  }
1571  return false;
1572 }
1573 /* End of JLai Modifications */
int getType() const
Real getCharge() const
void getBond(int num, int *atomi, int *atomj, int *bondtype) const
GenericMol(const char *theName)
int getNumBonds() const
int addPairGaussType2(int typea, int typeb, Real gA, Real gMu1, Real gSigma1)
GenericBond(int i, int j, int theType)
int getParams(const char *typea, const char *typeb, Real *c6, Real *c12, Real *c6pair, Real *c12pair) const
void getPairLJArrays2(int *indexA, int *indexB, Real *pairC6, Real *pairC12)
static bool GroGaussCompare(GroGaussPair, GroGaussPair)
int numGaussPair
int getNumDihedrals() const
int getNumDihedrals() const
const GenericBond * getBond(int n) const
#define BONDTYPES
const char * getAtomName() const
void add14Type(const char *typea, const char *typeb, Real c6pair, Real c12pair)
int getAtoml() const
int getResNum() const
bool bool_negative_number_warning_flag
void getVDWParams(int typea, int typeb, Real *c6, Real *c12, Real *c6pair, Real *c7) const
void getType(int num, char *type) const
const BigReal A
void addType(const char *type, Real m, Real q, Real c6, Real c12)
#define MOLECULES
int getParams(const char *type, Real *m, Real *q, Real *c6, Real *c12) const
ResizeArray< int > exclusions_atom_j
float Real
Definition: common.h:109
int getNumRes() const
#define ATOMS
int getNumAngles() const
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
#define iout
Definition: InfoStream.h:87
const char * getType() const
int getNumAtoms() const
#define NONBOND
#define DEFAULTS
#define ATOMTYPES
int getNumExclusions() const
const GenericAngle * getAngle(int n) const
void addType(const char *typea, const char *typeb, const char *typec, Real th0, Real kth, int funct)
ResizeArray< int > exclusions_atom_i
int getNum() const
#define BONDS
void getBondParams(int num, Real *b0, Real *kB, int *funct) const
int getIndex(Real b0, Real kB, int funct)
int getNumLJPair() const
#define PI
Definition: common.h:83
#define MOLECULETYPE
void getDihedral(int num, int *atomi, int *atomj, int *atomk, int *atoml, int *type) const
#define PAIRS
int getIndex(Real th0, Real kth, int funct)
#define EXCLUSIONS
int getNumAngles() const
void addType(const char *typea, const char *typeb, const Real *c, int mult, int funct)
int getType() const
int getNumGaussPair() const
iterator end(void)
Definition: ResizeArray.h:37
int getAtomk() const
#define ANGLES
#define SYSTEM
const char * getResType() const
int getType() const
void getParams(int num, Real *b0, Real *kB, int *funct) const
void getParams(int num, Real *c, int *mult, int *funct) const
GromacsTopFile(char *filename)
void getExclusions(int *, int *) const
int addPairLJType2(int typea, int typeb, Real c6, Real c12)
void addAngle(int atomi, int atomj, int atomk, int type)
void getParams(int num, Real *th0, Real *kth, int *funct) const
#define ANGLETYPES
int getNumBonds() const
void NAMD_die(const char *err_msg)
Definition: common.C:83
const GenericAtom * getAtom(int n) const
#define NAMESIZE
Definition: GromacsTopFile.h:6
const GenericMol * getMol() const
void addType(const char *typea, const char *typeb, Real c6, Real c12)
void addType(const char *typea, const char *typeb, Real b0, Real kB, int funct)
#define JOULES_PER_CALORIE
int add(const Elem &elem)
Definition: ResizeArray.h:97
int getNumAngles() const
int numLJPair
int getTypeNum() const
BlockRadixSort::TempStorage sort
void getPairLJArrays2(int *indexA, int *indexB, Real *pairC6, Real *pairC12)
const GenericDihedral * getDihedral(int n) const
#define DIHEDRALTYPES
#define LINESIZE
int getIndex(const Real *c, int mult, int funct)
int numExclusion
void getAngle(int num, int *atomi, int *atomj, int *atomk, int *angletype) const
int getNumAtoms() const
int getAtomj() const
void getDihedralParams(int num, Real *c, int *mult, int *funct) const
static bool GroLJCompare(GroLJPair, GroLJPair)
GenericAngle(int i, int j, int k, int theType)
int getAtomj() const
const BigReal B
int getAtomi() const
void getAtomParams(int num, char *type) const
#define LONGNAMESIZE
Definition: GromacsTopFile.h:7
k< npairi;++k){TABENERGY(const int numtypes=simParams->tableNumTypes;const float table_spacing=simParams->tableSpacing;const int npertype=(int)(namdnearbyint(simParams->tableMaxDist/simParams->tableSpacing)+1);) int table_i=(r2iilist[2 *k] >> 14)+r2_delta_expc;const int j=pairlisti[k];#define p_j BigReal diffa=r2list[k]-r2_table[table_i];#define table_four_i TABENERGY(register const int tabtype=-1-(lj_pars->A< 0?lj_pars->A:0);) BigReal kqq=kq_i *p_j-> charge
int getNumRes() const
Real getMass() const
void getPairGaussArrays2(int *indexA, int *indexB, Real *gaussA, Real *gaussMu1, Real *gaussSigma1, Real *gaussMu2, Real *gaussSigma2, Real *gaussRepulsive)
int getNumAtoms() const
int getNumDihedrals() const
void getPairGaussArrays2(int *indexA, int *indexB, Real *gaussA, Real *gaussMu1, Real *gaussSigma1, Real *gaussMu2, Real *gaussSigma2, Real *gaussRepulsive)
int size(void) const
Definition: ResizeArray.h:127
#define UNKNOWN
int getAtomi() const
int getAtomk() const
int getAtomi() const
void addBond(int atomi, int atomj, int type)
#define DIHEDRALS
int getAtomj() const
void addDihedral(int atomi, int atomj, int atomk, int atoml, int type)
int getNumPair() const
#define ANGSTROMS_PER_NM
GenericAtom(const char *theType, int theTypeNum, int theResNum, const char *theResType, const char *theAtomName, Real theCharge, Real theMass)
GenericDihedral(int i, int j, int k, int l, int theType)
void addAtom(const char *theType, int theTypeNum, int theResNum, const char *theResType, const char *theAtomName, Real theCharge, Real theMass)
void getAtom(int num, int *residue_number, char *residue_name, char *atom_name, char *atom_type, int *atom_typenum, Real *charge, Real *mass) const
void getAngleParams(int num, Real *th0, Real *kth, int *funct) const
iterator begin(void)
Definition: ResizeArray.h:36
MolInst(const GenericMol *theMol, int theNum)
int getNumBonds() const