NAMD
Public Member Functions | List of all members
GromacsTopFile Class Reference

#include <GromacsTopFile.h>

Public Member Functions

 GromacsTopFile (char *filename)
 
char * getSystemName () const
 
int getNumAtoms () const
 
int getNumBonds () const
 
int getNumAngles () const
 
int getNumDihedrals () const
 
int getNumAtomParams () const
 
int getNumBondParams () const
 
int getNumAngleParams () const
 
int getNumDihedralParams () const
 
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 getAtomParams (int num, char *type) const
 
int getNumPair () const
 
int getNumLJPair () const
 
int getNumGaussPair () const
 
int getNumExclusions () const
 
void getExclusions (int *, int *) const
 
void getBond (int num, int *atomi, int *atomj, int *bondtype) const
 
void getBondParams (int num, Real *b0, Real *kB, int *funct) const
 
void getAngle (int num, int *atomi, int *atomj, int *atomk, int *angletype) const
 
void getAngleParams (int num, Real *th0, Real *kth, int *funct) const
 
void getDihedral (int num, int *atomi, int *atomj, int *atomk, int *atoml, int *type) const
 
void getDihedralParams (int num, Real *c, int *mult, int *funct) const
 
void getPairLJArrays2 (int *indexA, int *indexB, Real *pairC6, Real *pairC12)
 
void getPairGaussArrays2 (int *indexA, int *indexB, Real *gaussA, Real *gaussMu1, Real *gaussSigma1, Real *gaussMu2, Real *gaussSigma2, Real *gaussRepulsive)
 
void getVDWParams (int typea, int typeb, Real *c6, Real *c12, Real *c6pair, Real *c7) const
 

Detailed Description

Definition at line 477 of file GromacsTopFile.h.

Constructor & Destructor Documentation

GromacsTopFile::GromacsTopFile ( char *  filename)

Definition at line 65 of file GromacsTopFile.C.

References ResizeArray< T >::add(), PairTable::addPairGaussType2(), PairTable::addPairLJType2(), AtomTable::addType(), BondTable::addType(), AngleTable::addType(), DihedralTable::addType(), VDWTable::addType(), ANGLES, ANGLETYPES, ANGSTROMS_PER_NM, ATOMS, ATOMTYPES, BONDS, BONDTYPES, bool_negative_number_warning_flag, charge, DEFAULTS, DIHEDRALS, DIHEDRALTYPES, EXCLUSIONS, BondTable::getIndex(), AngleTable::getIndex(), DihedralTable::getIndex(), AtomTable::getParams(), BondTable::getParams(), AngleTable::getParams(), DihedralTable::getParams(), iout, iWARN(), JOULES_PER_CALORIE, LINESIZE, LONGNAMESIZE, MOLECULES, MOLECULETYPE, NAMD_die(), NAMESIZE, NONBOND, numExclusion, PAIRS, PI, ResizeArray< T >::size(), SYSTEM, and UNKNOWN.

65  {
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 }
int addPairGaussType2(int typea, int typeb, Real gA, Real gMu1, Real gSigma1)
#define BONDTYPES
bool bool_negative_number_warning_flag
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:107
#define ATOMS
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
#define iout
Definition: InfoStream.h:87
#define NONBOND
#define DEFAULTS
#define ATOMTYPES
void addType(const char *typea, const char *typeb, const char *typec, Real th0, Real kth, int funct)
ResizeArray< int > exclusions_atom_i
#define BONDS
int getIndex(Real b0, Real kB, int funct)
#define PI
Definition: common.h:81
#define MOLECULETYPE
#define PAIRS
int getIndex(Real th0, Real kth, int funct)
#define EXCLUSIONS
void addType(const char *typea, const char *typeb, const Real *c, int mult, int funct)
#define ANGLES
#define SYSTEM
void getParams(int num, Real *b0, Real *kB, int *funct) const
void getParams(int num, Real *c, int *mult, int *funct) const
int addPairLJType2(int typea, int typeb, Real c6, Real c12)
void getParams(int num, Real *th0, Real *kth, int *funct) const
#define ANGLETYPES
void NAMD_die(const char *err_msg)
Definition: common.C:83
#define NAMESIZE
Definition: GromacsTopFile.h:6
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
#define DIHEDRALTYPES
#define LINESIZE
int getIndex(const Real *c, int mult, int funct)
int numExclusion
#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 size(void) const
Definition: ResizeArray.h:127
#define UNKNOWN
#define DIHEDRALS
#define ANGSTROMS_PER_NM

Member Function Documentation

void GromacsTopFile::getAngle ( int  num,
int *  atomi,
int *  atomj,
int *  atomk,
int *  angletype 
) const

Definition at line 727 of file GromacsTopFile.C.

References GenericAngle::getAtomi(), GenericAngle::getAtomj(), GenericAngle::getAtomk(), GenericAngle::getType(), and ResizeArray< T >::size().

728  {
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 }
int getAtomk() const
int getType() const
int getAtomi() const
int size(void) const
Definition: ResizeArray.h:127
int getAtomj() const
void GromacsTopFile::getAngleParams ( int  num,
Real th0,
Real kth,
int *  funct 
) const

Definition at line 1084 of file GromacsTopFile.C.

References AngleTable::getParams().

1085  {
1086  angleTable.getParams(num,th0,kth,funct);
1087 }
void getParams(int num, Real *th0, Real *kth, int *funct) const
void GromacsTopFile::getAtom ( int  num,
int *  residue_number,
char *  residue_name,
char *  atom_name,
char *  atom_type,
int *  atom_typenum,
Real charge,
Real mass 
) const

Definition at line 809 of file GromacsTopFile.C.

References GenericAtom::getAtomName(), GenericAtom::getCharge(), GenericAtom::getMass(), GenericAtom::getResNum(), GenericAtom::getResType(), GenericAtom::getType(), GenericAtom::getTypeNum(), and ResizeArray< T >::size().

Referenced by PDB::PDB().

813  {
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 }
Real getCharge() const
const char * getAtomName() const
int getResNum() const
const char * getType() const
const char * getResType() const
int getTypeNum() const
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
Real getMass() const
int size(void) const
Definition: ResizeArray.h:127
void GromacsTopFile::getAtomParams ( int  num,
char *  type 
) const
inline

Definition at line 542 of file GromacsTopFile.h.

References AtomTable::getType().

Referenced by getVDWParams().

542  {
543  atomTable.getType(num,type);
544  }
void getType(int num, char *type) const
void GromacsTopFile::getBond ( int  num,
int *  atomi,
int *  atomj,
int *  bondtype 
) const

Definition at line 693 of file GromacsTopFile.C.

References GenericBond::getAtomi(), GenericBond::getAtomj(), GenericBond::getType(), and ResizeArray< T >::size().

693  {
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 }
int getType() const
int getAtomj() const
int size(void) const
Definition: ResizeArray.h:127
int getAtomi() const
void GromacsTopFile::getBondParams ( int  num,
Real b0,
Real kB,
int *  funct 
) const

Definition at line 1077 of file GromacsTopFile.C.

References BondTable::getParams().

1078  {
1079  bondTable.getParams(num,b0,kB,funct);
1080 }
void getParams(int num, Real *b0, Real *kB, int *funct) const
void GromacsTopFile::getDihedral ( int  num,
int *  atomi,
int *  atomj,
int *  atomk,
int *  atoml,
int *  type 
) const

Definition at line 762 of file GromacsTopFile.C.

References GenericDihedral::getAtomi(), GenericDihedral::getAtomj(), GenericDihedral::getAtomk(), GenericDihedral::getAtoml(), GenericDihedral::getType(), and ResizeArray< T >::size().

763  {
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 }
int getType() const
int getAtoml() const
int getAtomj() const
int size(void) const
Definition: ResizeArray.h:127
int getAtomk() const
int getAtomi() const
void GromacsTopFile::getDihedralParams ( int  num,
Real c,
int *  mult,
int *  funct 
) const

Definition at line 1089 of file GromacsTopFile.C.

References DihedralTable::getParams().

1090  {
1091  dihedralTable.getParams(num,c,mult,funct);
1092 }
void getParams(int num, Real *c, int *mult, int *funct) const
void GromacsTopFile::getExclusions ( int *  atomi,
int *  atomj 
) const

Definition at line 680 of file GromacsTopFile.C.

References ResizeArray< T >::size().

680  {
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 }
ResizeArray< int > exclusions_atom_j
ResizeArray< int > exclusions_atom_i
int size(void) const
Definition: ResizeArray.h:127
int GromacsTopFile::getNumAngleParams ( ) const
inline

Definition at line 520 of file GromacsTopFile.h.

References AngleTable::size().

520 { return angleTable.size(); }
int size() const
int GromacsTopFile::getNumAngles ( ) const

Definition at line 639 of file GromacsTopFile.C.

References ResizeArray< T >::size().

639  {
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 }
int size(void) const
Definition: ResizeArray.h:127
int GromacsTopFile::getNumAtomParams ( ) const
inline

Definition at line 518 of file GromacsTopFile.h.

References AtomTable::size().

518 { return atomTable.size(); }
int size() const
int GromacsTopFile::getNumAtoms ( ) const

Definition at line 796 of file GromacsTopFile.C.

References ResizeArray< T >::size().

Referenced by PDB::PDB().

796  {
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 }
int size(void) const
Definition: ResizeArray.h:127
int GromacsTopFile::getNumBondParams ( ) const
inline

Definition at line 519 of file GromacsTopFile.h.

References BondTable::size().

519 { return bondTable.size(); }
int size() const
int GromacsTopFile::getNumBonds ( ) const

Definition at line 629 of file GromacsTopFile.C.

References ResizeArray< T >::size().

629  {
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 }
int size(void) const
Definition: ResizeArray.h:127
int GromacsTopFile::getNumDihedralParams ( ) const
inline

Definition at line 521 of file GromacsTopFile.h.

References DihedralTable::size().

521 { return dihedralTable.size(); }
int size() const
int GromacsTopFile::getNumDihedrals ( ) const

Definition at line 649 of file GromacsTopFile.C.

References ResizeArray< T >::size().

649  {
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 }
int size(void) const
Definition: ResizeArray.h:127
int GromacsTopFile::getNumExclusions ( ) const

Definition at line 675 of file GromacsTopFile.C.

References numExclusion.

675  {
676  return numExclusion;
677 }
int numExclusion
int GromacsTopFile::getNumGaussPair ( ) const

Definition at line 670 of file GromacsTopFile.C.

References numGaussPair.

670  {
671  return numGaussPair;
672 }
int numGaussPair
int GromacsTopFile::getNumLJPair ( ) const

Definition at line 666 of file GromacsTopFile.C.

References numLJPair.

666  {
667  return numLJPair;
668 }
int numLJPair
int GromacsTopFile::getNumPair ( ) const

Definition at line 660 of file GromacsTopFile.C.

References numGaussPair, and numLJPair.

660  {
661  int numPair = 0;
662  numPair = numLJPair + numGaussPair;
663  return numPair;
664 }
int numGaussPair
int numLJPair
void GromacsTopFile::getPairGaussArrays2 ( int *  indexA,
int *  indexB,
Real gaussA,
Real gaussMu1,
Real gaussSigma1,
Real gaussMu2,
Real gaussSigma2,
Real gaussRepulsive 
)

Definition at line 1098 of file GromacsTopFile.C.

References PairTable::getPairGaussArrays2().

1100  {
1101  pairTable.getPairGaussArrays2(indexA, indexB, gaussA, gaussMu1, gaussSigma1, gaussMu2, gaussSigma2, gaussRepulsive);
1102 }
void getPairGaussArrays2(int *indexA, int *indexB, Real *gaussA, Real *gaussMu1, Real *gaussSigma1, Real *gaussMu2, Real *gaussSigma2, Real *gaussRepulsive)
void GromacsTopFile::getPairLJArrays2 ( int *  indexA,
int *  indexB,
Real pairC6,
Real pairC12 
)

Definition at line 1094 of file GromacsTopFile.C.

References PairTable::getPairLJArrays2().

1094  {
1095  pairTable.getPairLJArrays2(indexA, indexB, pairC6, pairC12);
1096 }
void getPairLJArrays2(int *indexA, int *indexB, Real *pairC6, Real *pairC12)
char* GromacsTopFile::getSystemName ( ) const
inline

Definition at line 509 of file GromacsTopFile.h.

509 { return systemName; }
void GromacsTopFile::getVDWParams ( int  typea,
int  typeb,
Real c6,
Real c12,
Real c6pair,
Real c7 
) const

Definition at line 1417 of file GromacsTopFile.C.

References getAtomParams(), AtomTable::getParams(), and VDWTable::getParams().

1418  {
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 }
int getParams(const char *typea, const char *typeb, Real *c6, Real *c12, Real *c6pair, Real *c12pair) const
int getParams(const char *type, Real *m, Real *q, Real *c6, Real *c12) const
float Real
Definition: common.h:107
void getAtomParams(int num, char *type) const

The documentation for this class was generated from the following files: