NAMD
GoMolecule.C
Go to the documentation of this file.
1 
7 #include <stdio.h>
8 #include <string.h>
9 #include <stdlib.h>
10 #include <ctype.h>
11 #include "ResizeArray.h"
12 #include "InfoStream.h"
13 #include "Molecule.h"
14 #include "strlib.h"
15 #include "MStream.h"
16 #include "Communicate.h"
17 #include "Node.h"
18 #include "ObjectArena.h"
19 #include "Parameters.h"
20 #include "PDB.h"
21 #include "SimParameters.h"
22 #include "Hydrogen.h"
23 #include "UniqueSetIter.h"
24 #include "charm++.h"
25 /* BEGIN gf */
26 #include "ComputeGridForce.h"
27 #include "GridForceGrid.h"
28 #include "MGridforceParams.h"
29 /* END gf */
30 
31 #define MIN_DEBUG_LEVEL 3
32 //#define DEBUGM
33 #include "Debug.h"
34 #include "CompressPsf.h"
35 #include "ParallelIOMgr.h"
36 #include <deque>
37 #include <algorithm>
38 #ifndef M_PI
39 #define M_PI 3.14159265358979323846
40 #endif
41 #ifndef CODE_REDUNDANT
42 #define CODE_REDUNDANT 0
43 #endif
44 
45 // Ported by JLai from NAMD 2.7
46 /************************************************************************/
47 /* */
48 /* FUNCTION goInit */
49 /* */
50 /* This function is only called from the Molecule constructor. */
51 /* It only builds Go specific code into the Molecule object */
52 /* */
53 /************************************************************************/
54 void Molecule::goInit() {
55  numGoAtoms=0;
56  energyNative=0;
57  energyNonnative=0;
58  atomChainTypes=NULL;
59  goSigmaIndices=NULL;
60  goSigmas=NULL;
61  goWithinCutoff=NULL;
62  goCoordinates=NULL;
63  goResids=NULL;
64  goPDB=NULL;
65 }
66 
67 /************************************************************************/
68 /* */
69 /* FUNCTION build_go_params */
70 /* */
71 /* INPUTS: */
72 /* fname - name of the parameter file to read */
73 /* */
74 /* This function reads in multiple Go parameter files */
75 /* from a StringList and exhaustively processes them. */
76 /* */
77 /************************************************************************/
78 // JE - read in a Go parameter file
79 void Molecule::build_go_params(StringList *g) {
80  iout << iINFO << "Building Go Parameters" << "\n" << endi;
81 #ifdef MEM_OPT_VERSION
82  NAMD_die("Go forces are not supported in memory-optimized builds.");
83 #else
84  build_lists_by_atom();
85 #endif
86  int iterator = 0;
87  do
88  {
89  iout << iINFO << "Reading Go File: " << iterator << "\n" << endi;
90  read_go_file(g->data);
91  g = g->next;
92  iterator++;
93  } while ( g != NULL && iterator < 100);
94 }
95 
96 // Ported by JLai from NAMD 2.7
97 /************************************************************************/
98 /* */
99 /* FUNCTION read_go_file */
100 /* */
101 /* INPUTS: */
102 /* fname - name of the parameter file to read */
103 /* */
104 /* This function reads in a Go parameter file and adds the */
105 /* parameters from this file to the current group of parameters. */
106 /* The basic logic of the routine is to first find out what type of */
107 /* parameter we have in the file. Then look at each line in turn */
108 /* and call the appropriate routine to add the parameters until we hit*/
109 /* a new type of parameter or EOF. */
110 /* */
111 /************************************************************************/
112 // JE - read in a Go parameter file
113 void Molecule::read_go_file(char *fname)
114 
115 {
116 
117  int i; // Counter
118  int j; // Counter
119  int par_type=0; // What type of parameter are we currently
120  // dealing with? (vide infra)
121  // JLai -- uncommented
122  int skipline; // skip this line?
123  int skipall = 0; // skip rest of file;
124  char buffer[512]; // Buffer to store each line of the file
125  char first_word[512]; // First word of the current line
126  int read_count = 0; // Count of input parameters on a given line
127  int chain1 = 0; // First chain type for pair interaction
128  int chain2 = 0; // Second chain type for pair interaction
129  int int1; // First parameter int
130  int int2; // Second parameter int
131  Real r1; // Parameter Real
132  char in2[512]; // Second parameter word
133  GoValue *goValue1 = NULL; // current GoValue for loading parameters
134  GoValue *goValue2 = NULL; // other current GoValue for loading parameters
135  Bool sameGoChain = FALSE; // whether the current GoValue is within a chain or between chains
136  int restrictionCount = 0; // number of Go restrictions set for a given chain pair
137  FILE *pfile; // File descriptor for the parameter file
138 
139  /* Check to make sure that we haven't previously been told */
140  /* that all the files were read */
141  /*if (AllFilesRead)
142  {
143  NAMD_die("Tried to read another parameter file after being told that all files were read!");
144  }*/
145 
146  /* Initialize go_indices */
147  for (i=0; i<MAX_GO_CHAINS+1; i++) {
148  go_indices[i] = -1;
149  }
150 
151  /* Initialize go_array */
152  for (i=0; i<MAX_GO_CHAINS*MAX_GO_CHAINS; i++) {
153  go_array[i].epsilon = 1.25;
154  go_array[i].exp_a = 12;
155  go_array[i].exp_b = 6;
156  go_array[i].exp_rep = 12;
157  go_array[i].sigmaRep = 2.25;
158  go_array[i].epsilonRep = 0.03;
159  go_array[i].cutoff = 4.0;
160  for (j=0; j<MAX_RESTRICTIONS; j++) {
161  go_array[i].restrictions[j] = -1;
162  }
163  }
164 
165  /* Try and open the file */
166  if ( (pfile = fopen(fname, "r")) == NULL)
167  {
168  char err_msg[256];
169 
170  sprintf(err_msg, "UNABLE TO OPEN GO PARAMETER FILE %s\n", fname);
171  NAMD_die(err_msg);
172  }
173 
174  /* Keep reading in lines until we hit the EOF */
175  while (NAMD_read_line(pfile, buffer) != -1)
176  {
177  /* Get the first word of the line */
178  NAMD_find_first_word(buffer, first_word);
179  skipline=0;
180 
181  /* First, screen out things that we ignore. */
182  /* blank lines, lines that start with '!' or '*', lines that */
183  /* start with "END". */
184  if (!NAMD_blank_string(buffer) &&
185  (strncmp(first_word, "!", 1) != 0) &&
186  (strncmp(first_word, "*", 1) != 0) &&
187  (strncasecmp(first_word, "END", 3) != 0))
188  {
189  if ( skipall ) {
190  iout << iWARN << "SKIPPING PART OF GO PARAMETER FILE AFTER RETURN STATEMENT\n" << endi;
191  break;
192  }
193  /* Now, determine the apropriate parameter type. */
194  if (strncasecmp(first_word, "chaintypes", 10)==0)
195  {
196  read_count=sscanf(buffer, "%s %d %d\n", first_word, &int1, &int2);
197  if (read_count != 3) {
198  char err_msg[512];
199  sprintf(err_msg, "UNKNOWN PARAMETER IN GO PARAMETER FILE %s\nLINE=*%s*\nread_count=%d, int1=%d, int2=%d", fname, buffer, read_count, int1, int2);
200  NAMD_die(err_msg);
201  }
202  chain1 = int1;
203  chain2 = int2;
204  if (chain1 < 1 || chain1 > MAX_GO_CHAINS ||
205  chain2 < 1 || chain2 > MAX_GO_CHAINS) {
206  char err_msg[512];
207  sprintf(err_msg, "GO PARAMETER FILE: CHAIN INDEX MUST BE [1-%d] %s\nLINE=*%s*\nread_count=%d, int1=%d, int2=%d", MAX_GO_CHAINS, fname, buffer, read_count, int1, int2);
208  NAMD_die(err_msg);
209  }
210  if (go_indices[chain1] == -1) {
211  go_indices[chain1] = NumGoChains;
212  NumGoChains++;
213  }
214  if (go_indices[chain2] == -1) {
215  go_indices[chain2] = NumGoChains;
216  NumGoChains++;
217  }
218  if (chain1 == chain2) {
219  sameGoChain = TRUE;
220  } else {
221  sameGoChain = FALSE;
222  }
223  //goValue = &go_array[(chain1 * MAX_GO_CHAINS) + chain2];
224  goValue1 = &go_array[(chain1*MAX_GO_CHAINS) + chain2];
225  goValue2 = &go_array[(chain2*MAX_GO_CHAINS) + chain1];
226 #if CODE_REDUNDANT
227  goValue1 = &go_array[(go_indices[chain1]*MAX_GO_CHAINS) + go_indices[chain2]];
228  goValue2 = &go_array[(go_indices[chain2]*MAX_GO_CHAINS) + go_indices[chain1]];
229 #endif
230  restrictionCount = 0; // restrictionCount applies to each chain pair separately
231  }
232  else if (strncasecmp(first_word, "epsilonRep", 10)==0)
233  {
234  read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
235  if (read_count != 2) {}
236  goValue1->epsilonRep = r1;
237  if (!sameGoChain) {
238  goValue2->epsilonRep = r1;
239  }
240  }
241  else if (strncasecmp(first_word, "epsilon", 7)==0)
242  {
243  // Read in epsilon
244  read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
245  if (read_count != 2) {}
246  goValue1->epsilon = r1;
247  if (!sameGoChain) {
248  goValue2->epsilon = r1;
249  }
250  }
251  else if (strncasecmp(first_word, "exp_a", 5)==0)
252  {
253  read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
254  if (read_count != 2) {}
255  goValue1->exp_a = int1;
256  if (!sameGoChain) {
257  goValue2->exp_a = int1;
258  }
259  }
260  else if (strncasecmp(first_word, "exp_b", 5)==0)
261  {
262  read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
263  if (read_count != 2) {}
264  goValue1->exp_b = int1;
265  if (!sameGoChain) {
266  goValue2->exp_b = int1;
267  }
268  }
269  else if (strncasecmp(first_word, "exp_rep", 5)==0)
270  {
271  read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
272  if (read_count != 2) {}
273  goValue1->exp_rep = int1;
274  if (!sameGoChain) {
275  goValue2->exp_rep = int1;
276  }
277  }
278  else if (strncasecmp(first_word, "exp_Rep", 5)==0)
279  {
280  read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
281  if (read_count != 2) {}
282  goValue1->exp_rep = int1;
283  if (!sameGoChain) {
284  goValue2->exp_rep = int1;
285  }
286  }
287  else if (strncasecmp(first_word, "sigmaRep", 8)==0)
288  {
289  read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
290  if (read_count != 2) {}
291  goValue1->sigmaRep = r1;
292  if (!sameGoChain) {
293  goValue2->sigmaRep = r1;
294  }
295  }
296  else if (strncasecmp(first_word, "cutoff", 6)==0)
297  {
298  read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
299  if (read_count != 2) {}
300  goValue1->cutoff = r1;
301  if (!sameGoChain) {
302  goValue2->cutoff = r1;
303  }
304  }
305  else if (strncasecmp(first_word, "restriction", 10)==0)
306  {
307  read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
308  if (read_count != 2) {}
309  if (int1 < 0) {
310  DebugM(3, "ERROR: residue restriction value must be nonnegative. int1=" << int1 << "\n");
311  }
312  if (!sameGoChain) {
313  //goValue2->restrictions[restrictionCount] = int1;
314  DebugM(3, "ERROR: residue restrictions should not be defined between two separate chains. chain1=" << chain1 << ", chain2=" << chain2 << "\n");
315  }
316  else {
317  goValue1->restrictions[restrictionCount] = int1;
318  }
319  restrictionCount++;
320  }
321  else if (strncasecmp(first_word, "return", 4)==0)
322  {
323  skipall=8;
324  skipline=1;
325  }
326  else // if (par_type == 0)
327  {
328  /* This is an unknown paramter. */
329  /* This is BAD */
330  char err_msg[512];
331 
332  sprintf(err_msg, "UNKNOWN PARAMETER IN GO PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
333  NAMD_die(err_msg);
334  }
335  }
336  else
337  {
338  skipline=1;
339  }
340  }
341 
342  /* Close the file */
343  fclose(pfile);
344 
345  return;
346 }
347 /* END OF FUNCTION read_go_file */
348 
349 #if CODE_REDUNDANT
350 /************************************************************************/
351 /* */
352 /* FUNCTION get_go_cutoff */
353 /* */
354 /* INPUTS: */
355 /* chain1 - first chain type */
356 /* chain2 - second chain type */
357 /* */
358 /* This function gets the Go cutoff for two chain types. The cutoff */
359 /* determines whether the Go force is attractive or repulsive. */
360 /* */
361 /************************************************************************/
362 // JE
363 Real Molecule::get_go_cutoff(int chain1, int chain2)
364 {
365  return go_array[MAX_GO_CHAINS*chain1 + chain2].cutoff;
366  #if CODE_REDUNDANT
367  return go_array[MAX_GO_CHAINS*go_indices[chain1] + go_indices[chain2]].cutoff;
368  #endif
369 }
370 /* END OF FUNCTION get_go_cutoff */
371 
372 
373 /************************************************************************/
374 /* */
375 /* FUNCTION get_go_epsilonRep */
376 /* */
377 /* INPUTS: */
378 /* chain1 - first chain type */
379 /* chain2 - second chain type */
380 /* */
381 /* This function gets the Go epsilonRep value for two chain types. */
382 /* epsilonRep is a factor in the repulsive Go force formula. */
383 /* */
384 /************************************************************************/
385 // JE
386 Real Molecule::get_go_epsilonRep(int chain1, int chain2)
387 {
388  return go_array[MAX_GO_CHAINS*chain1 + chain2].epsilonRep;
389  #if CODE_REDUNDANT
390  return go_array[MAX_GO_CHAINS*go_indices[chain1] + go_indices[chain2]].epsilonRep;
391  #endif
392 }
393 /* END OF FUNCTION get_go_epsilonRep */
394 
395 
396 /************************************************************************/
397 /* */
398 /* FUNCTION get_go_sigmaRep */
399 /* */
400 /* INPUTS: */
401 /* chain1 - first chain type */
402 /* chain2 - second chain type */
403 /* */
404 /* This function gets the Go sigmaRep value for two chain types. */
405 /* sigmaRep is a factor in the repulsive Go force formula. */
406 /* */
407 /************************************************************************/
408 // JE
409 Real Molecule::get_go_sigmaRep(int chain1, int chain2)
410 {
411  return go_array[MAX_GO_CHAINS*chain1 + chain2].sigmaRep;
412 }
413 /* END OF FUNCTION get_go_sigmaRep */
414 
415 
416 /************************************************************************/
417 /* */
418 /* FUNCTION get_go_epsilon */
419 /* */
420 /* INPUTS: */
421 /* chain1 - first chain type */
422 /* chain2 - second chain type */
423 /* */
424 /* This function gets the Go epsilon value for two chain types. */
425 /* epsilon is a factor in the attractive Go force formula. */
426 /* */
427 /************************************************************************/
428 // JE
429 Real Molecule::get_go_epsilon(int chain1, int chain2)
430 {
431  return go_array[MAX_GO_CHAINS*chain1 + chain2].epsilon;
432  #if CODE_REDUNDANT
433  return go_array[MAX_GO_CHAINS*go_indices[chain1] + go_indices[chain2]].epsilon;
434  #endif
435 }
436 /* END OF FUNCTION get_go_epsilon */
437 
438 
439 /************************************************************************/
440 /* */
441 /* FUNCTION get_go_exp_a */
442 /* */
443 /* INPUTS: */
444 /* chain1 - first chain type */
445 /* chain2 - second chain type */
446 /* */
447 /* This function gets the Go exp_a value for two chain types. */
448 /* exp_a is an exponent in the repulsive term of the attractive Go */
449 /* force formula. */
450 /* */
451 /************************************************************************/
452 // JE
453 int Molecule::get_go_exp_a(int chain1, int chain2)
454 {
455  return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_a;
456  #if CODE_REDUNDANT
457  return go_array[MAX_GO_CHAINS*go_indices[chain1] + go_indices[chain2]].exp_a;
458  #endif
459 }
460 /* END OF FUNCTION get_go_exp_a */
461 
462 
463 /************************************************************************/
464 /* */
465 /* FUNCTION get_go_exp_b */
466 /* */
467 /* INPUTS: */
468 /* chain1 - first chain type */
469 /* chain2 - second chain type */
470 /* */
471 /* This function gets the Go exp_b value for two chain types. */
472 /* exp_b is an exponent in the attractive term of the attractive Go */
473 /* force formula. */
474 /* */
475 /************************************************************************/
476 // JE
477 int Molecule::get_go_exp_b(int chain1, int chain2)
478 {
479  return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_b;
480  #if CODE_REDUNDANT
481  return go_array[MAX_GO_CHAINS*go_indices[chain1] + go_indices[chain2]].exp_b;
482  #endif
483 }
484 /* END OF FUNCTION get_go_exp_b */
485 
486 
487 /************************************************************************/
488 /* */
489 /* FUNCTION get_go_exp_rep */
490 /* */
491 /* INPUTS: */
492 /* chain1 - first chain type */
493 /* chain2 - second chain type */
494 /* */
495 /* This function gets the Go exp_rep value for two chain types. */
496 /* exp_b is an exponent in the attractive term of the attractive Go */
497 /* force formula. */
498 /* */
499 /************************************************************************/
500 // JE
501 int Molecule::get_go_exp_rep(int chain1, int chain2)
502 {
503  return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_rep;
504  #if CODE_REDUNDANT
505  return go_array[MAX_GO_CHAINS*go_indices[chain1] + go_indices[chain2]].exp_rep;
506  #endif
507 }
508 /* END OF FUNCTION get_go_exp_rep */
509 #endif
510 
511 /************************************************************************/
512 /* */
513 /* FUNCTION go_restricted */
514 /* */
515 /* INPUTS: */
516 /* chain1 - first chain type */
517 /* chain2 - second chain type */
518 /* rDiff - residue ID difference to check */
519 /* */
520 /* This function checks whether residues with IDs rDiff apart are */
521 /* restricted from Go force calculation. */
522 /* */
523 /************************************************************************/
524 // JE
525 Bool Molecule::go_restricted(int chain1, int chain2, int rDiff)
526 {
527  int i; // Loop counter
528  for(i=0; i<MAX_RESTRICTIONS; i++) {
529  if (go_array[(MAX_GO_CHAINS*chain1) + chain2].restrictions[i] == rDiff) {
530  return TRUE;
531  } else if (go_array[(MAX_GO_CHAINS*chain1) + chain2].restrictions[i] == -1) {
532  return FALSE;
533  }
534  }
535  return FALSE;
536 }
537 /* END OF FUNCTION go_restricted */
538 
539 // Original by JE
540 /************************************************************************/
541 /* */
542 /* FUNCTION print_go_params */
543 /* */
544 /* This is a debugging routine used to print out all the Go */
545 /* parameters */
546 /* */
547 /************************************************************************/
548 void Molecule::print_go_params()
549 {
550  int i;
551  int j;
552  int index;
553 
554  DebugM(3,NumGoChains << " Go PARAMETERS 3\n" \
555  << "*****************************************" << std::endl);
556 
557  for (i=0; i<NumGoChains; i++) {
558  for (j=0; j<NumGoChains; j++) {
559  index = (i * MAX_GO_CHAINS) + j;
560  // Real epsilon; // Epsilon
561  // Real exp_a; // First exponent for attractive L-J term
562  // Real exp_b; // Second exponent for attractive L-J term
563  // Real sigmaRep; // Sigma for repulsive term
564  // Real epsilonRep; // Epsilon for replusive term
565  DebugM(3,"Go index=(" << i << "," << j << ") epsilon=" << go_array[index].epsilon \
566  << " exp_a=" << go_array[index].exp_a << " exp_b=" << go_array[index].exp_b \
567  << " exp_rep=" << go_array[index].exp_rep << " sigmaRep=" \
568  << go_array[index].sigmaRep << " epsilonRep=" << go_array[index].epsilonRep \
569  << " cutoff=" << go_array[index].cutoff << std::endl);
570  }
571  }
572 
573 }
574 // End of port -- JLai
575 
576 #ifndef MEM_OPT_VERSION
577 void Molecule::build_go_sigmas(StringList *goCoordFile,
578  char *cwd)
579 {
580  DebugM(3,"->build_go_sigmas" << std::endl);
581  PDB *goPDB; // Pointer to PDB object to use
582  int bcol = 4; // Column that data is in
583  int32 chainType = 0; // b value from PDB file
584  int i; // Loop counter
585  int j; // Loop counter
586  int resid1; // Residue ID for first atom
587  int resid2; // Residue ID for second atom
588  int residDiff; // Difference between resid1 and resid2
589  Real sigma; // Sigma calculated for a Go atom pair
590  Real atomAtomDist; // Distance between two atoms
591  Real exp_a; // First exponent in L-J like formula
592  Real exp_b; // Second exponent in L-J like formula
593  char filename[NAMD_FILENAME_BUFFER_SIZE]; // Filename
594 
595  //JLai
596  BigReal nativeEnergy, *native;
597  BigReal nonnativeEnergy, *nonnative;
598  nativeEnergy = 0;
599  nonnativeEnergy = 0;
600  native = &nativeEnergy;
601  nonnative = &nonnativeEnergy;
602 
603  long nativeContacts = 0;
604  long nonnativeContacts = 0;
605 
606  // Get the PDB object that contains the Go coordinates. If
607  // the user gave another file name, use it. Otherwise, just use
608  // the PDB file that has the initial coordinates.
609  if (goCoordFile == NULL)
610  {
611  //goPDB = initial_pdb;
612  NAMD_die("Error: goCoordFile is NULL - build_go_sigmas");
613  }
614  else
615  {
616  if (goCoordFile->next != NULL)
617  {
618  NAMD_die("Multiple definitions of Go atoms PDB file in configuration file");
619  }
620 
621  if ( (cwd == NULL) || (goCoordFile->data[0] == '/') )
622  {
623  strcpy(filename, goCoordFile->data);
624  }
625  else
626  {
627  strcpy(filename, cwd);
628  strcat(filename, goCoordFile->data);
629  }
630 
631  goPDB = new PDB(filename);
632  if ( goPDB == NULL )
633  {
634  NAMD_die("Memory allocation failed in Molecule::build_go_sigmas");
635  }
636 
637  if (goPDB->num_atoms() != numAtoms)
638  {
639  NAMD_die("Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
640  }
641  }
642  // Allocate the array to hold the chain types
643  atomChainTypes = new int32[numAtoms];
644  // Allocate the array to hold Go atom indices into the sigma array
645  goSigmaIndices = new int32[numAtoms];
646 
647  if (atomChainTypes == NULL) {
648  NAMD_die("memory allocation failed in Molecule::build_go_sigmas");
649  }
650 
651  numGoAtoms = 0;
652 
653  // Loop through all the atoms and get the Go chain types
654  for (i=0; i<numAtoms; i++) {
655  // Get the chainType from the occupancy field
656  chainType = (int32)(goPDB->atom(i))->occupancy();
657  // Assign the chainType value
658  if ( chainType != 0 ) {
659  //DebugM(3,"build_go_sigmas - atom:" << i << ", chainType:" << chainType << std::endl);
660  atomChainTypes[i] = chainType;
661  goSigmaIndices[i] = numGoAtoms;
662  numGoAtoms++;
663  }
664  else {
665  atomChainTypes[i] = 0;
666  goSigmaIndices[i] = -1;
667  }
668  //printf("CT: %d %d %d %d\n",i,numGoAtoms,atomChainTypes[i],goSigmaIndices[i]);
669  }
670 
671  // Allocate the array to hold the sigma values for each Go atom pair
672  goSigmas = new Real[numGoAtoms*numGoAtoms];
673  goWithinCutoff = new bool[numGoAtoms*numGoAtoms];
674  for (i=0; i<numGoAtoms; i++) {
675  for (j=0; j<numGoAtoms; j++) {
676  goSigmas[i*numGoAtoms + j] = -1.0;
677  goWithinCutoff[i*numGoAtoms + j] = false;
678  }
679  }
680  // Loop through all atom pairs and calculate sigmas
681  DebugM(3," numAtoms=" << numAtoms << std::endl);
682  for (i=0; i<numAtoms; i++) {
683  //DebugM(3," i=" << i << std::endl);
684  resid1 = (goPDB->atom(i))->residueseq();
685  //DebugM(3," resid1=" << resid1 << std::endl);
686  //if ( goSigmaIndices[i] != -1) {
687  // goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[i]] = 0.0;
688  //}
689  for (j=i+1; j<numAtoms; j++) {
690  //DebugM(3," j=" << j << std::endl);
691  resid2 = (goPDB->atom(j))->residueseq();
692  //printf("GSIi %d %d %d\n",i,numAtoms,goSigmaIndices[i]);
693  //printf("SAN CHECK: %d\n",goSigmaIndices[37]);
694  //printf("GSIj %d %d %d\n",j,numAtoms,goSigmaIndices[j]);
695  //printf("ATOMS_1to4 %d\n",!atoms_1to4(i,j));
696  //DebugM(3," resid2=" << resid2 << std::endl);
697  // if goSigmaIndices aren't defined, don't set anything in goSigmas
698  if ( goSigmaIndices[i] != -1 && goSigmaIndices[j] != -1 && !atoms_1to4(i,j) ) {
699  //printf("TAKING DIFFERENCE\n");
700  residDiff = resid2 - resid1;
701  //printf("RESIDDIFF %d\n",residDiff);
702  if (residDiff < 0) residDiff = -residDiff;
703  //printf("RESIDDIFF2 %d\n",residDiff);
704  // if this is a Go atom pair that is not restricted
705  // calculate sigma
706  // sigmas are initially set to -1.0 if the atom pair fails go_restricted
707  //printf("CHECKING LOOPING\n");
708  if ( atomChainTypes[i] && atomChainTypes[j] &&
709  !(this->go_restricted(atomChainTypes[i],atomChainTypes[j],residDiff)) ) {
710  atomAtomDist = sqrt(pow((goPDB->atom(i))->xcoor() - (goPDB->atom(j))->xcoor(), 2.0) +
711  pow((goPDB->atom(i))->ycoor() - (goPDB->atom(j))->ycoor(), 2.0) +
712  pow((goPDB->atom(i))->zcoor() - (goPDB->atom(j))->zcoor(), 2.0));
713  if ( atomAtomDist <= this->get_go_cutoff(atomChainTypes[i],atomChainTypes[j]) ) {
714  exp_a = this->get_go_exp_a(atomChainTypes[i],atomChainTypes[j]);
715  exp_b = this->get_go_exp_b(atomChainTypes[i],atomChainTypes[j]);
716  sigma = pow(static_cast<double>(exp_b/exp_a),(1.0/(exp_a-exp_b))) * atomAtomDist;
717  goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = sigma;
718  goSigmas[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = sigma;
719  goWithinCutoff[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = true;
720  goWithinCutoff[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = true;
721  nativeContacts++;
722  } else {
723  goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = 0.0;
724  goSigmas[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = 0.0;
725  nonnativeContacts++;
726  }
727  } else {
728  goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = -1.0;
729  goSigmas[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = -1.0;
730  }
731  }
732  }
733  }
734 
735  iout << iINFO << "Number of UNIQUE native contacts: " << nativeContacts << "\n" << endi;
736  iout << iINFO << "Number of UNIQUE nonnative contacts: " << nonnativeContacts << "\n" << endi;
737 
738  // If we had to create a PDB object, delete it now
739  if (goCoordFile != NULL) {
740  delete goPDB;
741  }
742 
743  return;
744 }
745 /* END OF FUNCTION build_go_sigmas */
746 
747 void Molecule::build_go_sigmas2(StringList *goCoordFile,
748  char *cwd)
749 {
750  DebugM(3,"->build_go_sigmas2" << std::endl);
751  PDB *goPDB; // Pointer to PDB object to use
752  int bcol = 4; // Column that data is in
753  int32 chainType = 0; // b value from PDB file
754  int32 residType = 0; // resid value from PDB file
755  int i; // Loop counter
756  int j; // Loop counter
757  int resid1; // Residue ID for first atom
758  int resid2; // Residue ID for second atom
759  int residDiff; // Difference between resid1 and resid2
760  Real sigma; // Sigma calculated for a Go atom pair
761  Real atomAtomDist; // Distance between two atoms
762  Real exp_a; // First exponent in L-J like formula
763  Real exp_b; // Second exponent in L-J like formula
764  char filename[NAMD_FILENAME_BUFFER_SIZE]; // Filename
765 
766  //JLai
767  int numLJPair = 0;
768  long nativeContacts = 0;
769  long nonnativeContacts = 0;
770 
771  // Get the PDB object that contains the Go coordinates. If
772  // the user gave another file name, use it. Otherwise, just use
773  // the PDB file that has the initial coordinates.
774  if (goCoordFile == NULL)
775  {
776  //goPDB = initial_pdb;
777  NAMD_die("Error: goCoordFile is NULL - build_go_sigmas2");
778  }
779  else
780  {
781  if (goCoordFile->next != NULL)
782  {
783  NAMD_die("Multiple definitions of Go atoms PDB file in configuration file");
784  }
785 
786  if ( (cwd == NULL) || (goCoordFile->data[0] == '/') )
787  {
788  strcpy(filename, goCoordFile->data);
789  }
790  else
791  {
792  strcpy(filename, cwd);
793  strcat(filename, goCoordFile->data);
794  }
795 
796  goPDB = new PDB(filename);
797  if ( goPDB == NULL )
798  {
799  NAMD_die("Memory allocation failed in Molecule::build_go_sigmas2");
800  }
801 
802  if (goPDB->num_atoms() != numAtoms)
803  {
804  NAMD_die("Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
805  }
806  }
807  // Allocate the array to hold the chain types
808  atomChainTypes = new int32[numAtoms];
809  // Allocate the array to hold Go atom indices into the sigma array
810  goSigmaIndices = new int32[numAtoms];
811  // Allocate the array to hold resid
812  goResidIndices = new int32[numAtoms];
813 
814  if (atomChainTypes == NULL) {
815  NAMD_die("memory allocation failed in Molecule::build_go_sigmas2");
816  }
817 
818  numGoAtoms = 0;
819 
820  // Loop through all the atoms and get the Go chain types
821  for (i=0; i<numAtoms; i++) {
822  // Get the chainType from the occupancy field
823  chainType = (int32)(goPDB->atom(i))->occupancy();
824  // Get the resid from the resid field
825  residType = (int32)(goPDB->atom(i))->residueseq();
826  // Assign the chainType value
827  if ( chainType != 0 ) {
828  //DebugM(3,"build_go_sigmas2 - atom:" << i << ", chainType:" << chainType << std::endl);
829  atomChainTypes[i] = chainType;
830  goSigmaIndices[i] = numGoAtoms;
831  goResidIndices[i] = residType;
832  numGoAtoms++;
833  }
834  else {
835  atomChainTypes[i] = 0;
836  goSigmaIndices[i] = -1;
837  goResidIndices[i] = -1;
838  }
839  }
840 
841  //Define:
842  ResizeArray<GoPair> tmpGoPair;
843 
844  // Loop through all atom pairs and calculate sigmas
845  // Store sigmas into sorted array
846  DebugM(3," numAtoms=" << numAtoms << std::endl);
847  for (i=0; i<numAtoms; i++) {
848  resid1 = (goPDB->atom(i))->residueseq();
849  for (j=i+1; j<numAtoms; j++) {
850  resid2 = (goPDB->atom(j))->residueseq();
851  if ( goSigmaIndices[i] != -1 && goSigmaIndices[j] != -1 && !atoms_1to4(i,j) ) {
852  residDiff = resid2 - resid1;
853  if (residDiff < 0) residDiff = -residDiff;
854  if ( atomChainTypes[i] && atomChainTypes[j] &&
855  !(this->go_restricted(atomChainTypes[i],atomChainTypes[j],residDiff)) ) {
856  atomAtomDist = sqrt(pow((goPDB->atom(i))->xcoor() - (goPDB->atom(j))->xcoor(), 2.0) +
857  pow((goPDB->atom(i))->ycoor() - (goPDB->atom(j))->ycoor(), 2.0) +
858  pow((goPDB->atom(i))->zcoor() - (goPDB->atom(j))->zcoor(), 2.0));
859  if ( atomAtomDist <= this->get_go_cutoff(atomChainTypes[i],atomChainTypes[j]) ) {
860  exp_a = this->get_go_exp_a(atomChainTypes[i],atomChainTypes[j]);
861  exp_b = this->get_go_exp_b(atomChainTypes[i],atomChainTypes[j]);
862  sigma = pow(static_cast<double>(exp_b/exp_a),(1.0/(exp_a-exp_b))) * atomAtomDist;
863  double tmpA = pow(sigma,exp_a);
864  double tmpB = pow(sigma,exp_b);
865  GoPair gp;
866  GoPair gp2;
867  gp.goIndxA = i;
868  gp.goIndxB = j;
869  gp.A = tmpA;
870  gp.B = tmpB;
871  tmpGoPair.add(gp);
872  gp2.goIndxA = j;
873  gp2.goIndxB = i;
874  gp2.A = tmpA;
875  gp2.B = tmpB;
876  tmpGoPair.add(gp2);
877  nativeContacts++;
878  } else {
879  nonnativeContacts++;
880  }
881  }
882  }
883  }
884  }
885 
886  iout << iINFO << "Number of UNIQUE native contacts: " << nativeContacts << "\n" << endi;
887  iout << iINFO << "Number of UNIQUE nonnative contacts: " << nonnativeContacts << "\n" << endi;
888 
889  // Copies the resizeArray into a block of continuous memory
890  std::sort(tmpGoPair.begin(),tmpGoPair.end(),goPairCompare);
891  goNumLJPair = 2*nativeContacts;
892  goIndxLJA = new int[goNumLJPair];
893  goIndxLJB = new int[goNumLJPair];
894  goSigmaPairA = new double[goNumLJPair];
895  goSigmaPairB = new double[goNumLJPair];
896  for(i=0; i< goNumLJPair; i++) {
897  goIndxLJA[i] = tmpGoPair[i].goIndxA;
898  goIndxLJB[i] = tmpGoPair[i].goIndxB;
899  goSigmaPairA[i] = tmpGoPair[i].A;
900  goSigmaPairB[i] = tmpGoPair[i].B;
901  }
902 
903  pointerToGoBeg = new int[numAtoms];
904  pointerToGoEnd = new int[numAtoms];
905  int oldIndex = -1;
906  for(i=0; i<numAtoms; i++) {
907  pointerToGoBeg[i] = -1;
908  pointerToGoEnd[i] = -2;
909  }
910  for(i=0; i < goNumLJPair; i++) {
911  if(pointerToGoBeg[goIndxLJA[i]] == -1) {
912  pointerToGoBeg[goIndxLJA[i]] = i;
913  oldIndex = goIndxLJA[i];
914  }
915  pointerToGoEnd[oldIndex] = i;
916  }
917 
918  // If we had to create a PDB object, delete it now
919  if (goCoordFile != NULL) {
920  delete goPDB;
921  }
922 
923  return;
924 }
925 /* END OF FUNCTION build_go_sigmas2 */
926 
927 bool Molecule::goPairCompare(GoPair first, GoPair second) {
928  if(first.goIndxA < second.goIndxA) {
929  return true;
930  } else if(first.goIndxA == second.goIndxA) {
931  return (first.goIndxB == second.goIndxB);
932  }
933  return false;
934 }
935 
936  /************************************************************************/
937  /* */
938  /* JE - FUNCTION build_go_arrays */
939  /* */
940  /* INPUTS: */
941  /* goCoordFile - Value of Go coordinate file from config file */
942  /* cwd - Current working directory */
943  /* */
944  /* This function builds arrays that support sigma calculation for L-J */
945  /* style Go forces. It takes the name of the PDB file. It then builds */
946  /* an array identifying atoms to which Go forces apply. */
947  /* */
948  /************************************************************************/
949 // JE
950 void Molecule::build_go_arrays(StringList *goCoordFile,
951  char *cwd)
952 {
953  DebugM(3,"->build_go_arrays" << std::endl);
954  //PDB *goPDB; // Pointer to PDB object to use
955  int bcol = 4; // Column that data is in
956  int32 chainType = 0; // b value from PDB file
957  int i; // Loop counter
958  int j; // Loop counter
959  BigReal atomAtomDist; // Distance between two atoms -- JLai put back
960  int resid1; // Residue ID for first atom
961  int resid2; // Residue ID for second atom
962  int residDiff; // Difference between resid1 and resid2
963  int goIndex; // Index into the goCoordinates array
964  int goIndx; // Index into the goCoordinates array
965  char filename[NAMD_FILENAME_BUFFER_SIZE]; // Filename
966 
967  //JLai
968  BigReal nativeEnergy, *native;
969  BigReal nonnativeEnergy, *nonnative;
970  nativeEnergy = 0;
971  nonnativeEnergy = 0;
972  native = &nativeEnergy;
973  nonnative = &nonnativeEnergy;
974 
975  long nativeContacts = 0;
976  long nonnativeContacts = 0;
977 
978  // Get the PDB object that contains the Go coordinates. If
979  // the user gave another file name, use it. Otherwise, just use
980  // the PDB file that has the initial coordinates.
981  if (goCoordFile == NULL)
982  {
983  //goPDB = initial_pdb;
984  NAMD_die("Error: goCoordFile is NULL - build_go_arrays");
985  }
986  else
987  {
988  if (goCoordFile->next != NULL)
989  {
990  NAMD_die("Multiple definitions of Go atoms PDB file in configuration file");
991  }
992 
993  if ( (cwd == NULL) || (goCoordFile->data[0] == '/') )
994  {
995  strcpy(filename, goCoordFile->data);
996  }
997  else
998  {
999  strcpy(filename, cwd);
1000  strcat(filename, goCoordFile->data);
1001  }
1002 
1003  goPDB = new PDB(filename);
1004  if ( goPDB == NULL )
1005  {
1006  NAMD_die("goPDB memory allocation failed in Molecule::build_go_arrays");
1007  }
1008 
1009  if (goPDB->num_atoms() != numAtoms)
1010  {
1011  NAMD_die("Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
1012  }
1013  }
1014 
1015  // Allocate the array to hold Go atom indices into the sigma array
1016  goSigmaIndices = new int32[numAtoms];
1017 
1018  if (goSigmaIndices == NULL) {
1019  NAMD_die("goSigmaIndices memory allocation failed in Molecule::build_go_arrays");
1020  }
1021 
1022  numGoAtoms = 0;
1023 
1024  // Loop through all the atoms and get the Go chain types
1025  for (i=0; i<numAtoms; i++) {
1026  chainType = (int32)(goPDB->atom(i))->occupancy();
1027  if ( chainType != 0 ) {
1028  DebugM(3,"build_go_arrays - atom:" << i << std::endl);
1029  goSigmaIndices[i] = numGoAtoms;
1030  numGoAtoms++;
1031  }
1032  else {
1033  goSigmaIndices[i] = -1;
1034  }
1035  }
1036 
1037  // Allocate the array to hold the sigma values for each Go atom pair
1049  // Allocate the array to hold the chain types
1050  atomChainTypes = new int32[numGoAtoms];
1051 
1052  if (atomChainTypes == NULL) {
1053  NAMD_die("atomChainTypes memory allocation failed in Molecule::build_go_arrays");
1054  }
1055 
1056  // Allocate the array to hold (x,y,z) coordinates for all of the Go atoms
1057  goCoordinates = new Real[numGoAtoms*3];
1058 
1059  if (goCoordinates == NULL) {
1060  NAMD_die("goCoordinates memory allocation failed in Molecule::build_go_arrays");
1061  }
1062 
1063  goResids = new int[numGoAtoms];
1064 
1065  // Allocate the array to hold PDB residu IDs for all of the Go atoms
1066  if (goResids == NULL) {
1067  NAMD_die("goResids memory allocation failed in Molecule::build_go_arrays");
1068  }
1069 
1070  for (i=0; i<numAtoms; i++) {
1071  goIndex = goSigmaIndices[i];
1072  if (goIndex != -1) {
1073  // Assign the chainType value!
1074  // Get the chainType from the occupancy field
1075  atomChainTypes[goIndex] = (int32)(goPDB->atom(i))->occupancy();
1076  goCoordinates[goIndex*3] = goPDB->atom(i)->xcoor();
1077  goCoordinates[goIndex*3 + 1] = goPDB->atom(i)->ycoor();
1078  goCoordinates[goIndex*3 + 2] = goPDB->atom(i)->zcoor();
1079  goResids[goIndex] = goPDB->atom(i)->residueseq();
1080  }
1081  }
1082  // JLai
1083  energyNative = 0;
1084  energyNonnative = 0;
1085  //printf("INIT ENERGY: (N) %f (NN) %f\n", energyNative, energyNonnative);
1086  for (i=0; i<numAtoms-1; i++) {
1087  goIndex = goSigmaIndices[i];
1088  if (goIndex != -1) {
1089  for (j=i+1; j<numAtoms;j++) {
1090  goIndx = goSigmaIndices[j];
1091  if (goIndx != -1) {
1092  resid1 = (goPDB->atom(i))->residueseq();
1093  resid2 = (goPDB->atom(j))->residueseq();
1094  residDiff = resid2 - resid1;
1095  if (residDiff < 0) residDiff = -residDiff;
1096  if (atomChainTypes[goIndex] && atomChainTypes[goIndx] &&
1097  !(this->go_restricted(atomChainTypes[goIndex],atomChainTypes[goIndx],residDiff)) &&
1098  !atoms_1to4(i,j)) {
1099  atomAtomDist = sqrt(pow((goPDB->atom(i))->xcoor() - (goPDB->atom(j))->xcoor(), 2.0) +
1100  pow((goPDB->atom(i))->ycoor() - (goPDB->atom(j))->ycoor(), 2.0) +
1101  pow((goPDB->atom(i))->zcoor() - (goPDB->atom(j))->zcoor(), 2.0));
1102  if (atomAtomDist <= this->get_go_cutoff(atomChainTypes[goIndex],atomChainTypes[goIndx]) ) {
1103  nativeContacts++;
1104  } else {
1105  nonnativeContacts++;
1106  }
1107  }
1108  }
1109  }
1110  }
1111  }
1112  iout << iINFO << "Number of UNIQUE native contacts: " << nativeContacts << "\n" << endi;
1113  iout << iINFO << "Number of UNIQUE nonnative contacts: " << nonnativeContacts << "\n" << endi;
1114 
1115  // If we had to create a PDB object, delete it now
1116  if (goCoordFile != NULL) {
1117  delete goPDB;
1118  }
1119 
1120  return;
1121 }
1122 /* END OF FUNCTION build_go_arrays */
1123 #endif // #ifndef MEM_OPT_VERSION
1124 
1125 /************************************************************************/
1126 /* */
1127 /* FUNCTION print_go_sigmas */
1128 /* */
1129 /* print_go_sigmas prints out goSigmas, the array of sigma parameters */
1130 /* used in the L-J type Go force calculations */
1131 /* */
1132 /************************************************************************/
1133 // JE
1134 void Molecule::print_go_sigmas()
1135 {
1136  int i; // Counter
1137  int j; // Counter
1138  Real sigma;
1139 
1140  DebugM(3,"GO SIGMA ARRAY\n" \
1141  << "***************************" << std::endl);
1142  DebugM(3, "numGoAtoms: " << numGoAtoms << std::endl);
1143 
1144  if (goSigmaIndices == NULL) {
1145  DebugM(3, "GO SIGMAS HAVE NOT BEEN BUILT" << std::endl);
1146  return;
1147  }
1148 
1149  for (i=0; i<numAtoms; i++) {
1150  for (j=0; j<numAtoms; j++) {
1151  if ( goSigmaIndices[i] != -1 && goSigmaIndices[j] != -1 ) {
1152  //DebugM(3, "i: " << i << ", j: " << j << std::endl);
1153  sigma = goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]];
1154  if (sigma > 0.0) {
1155  DebugM(3, "(" << i << "," << j << ") - +" << sigma << " ");
1156  }
1157  else {
1158  DebugM(3, "(" << i << "," << j << ") - " << sigma << " ");
1159  }
1160  } else {
1161  //DebugM(3, "0 ");
1162  }
1163  }
1164  if ( goSigmaIndices[i] != -1 ) {
1165  DebugM(3, "-----------" << std::endl);
1166  }
1167  }
1168  return;
1169 }
1170 /* END OF FUNCTION print_go_sigmas */
1171 
1172 // JLai
1173 BigReal Molecule::get_gro_force2(BigReal x,
1174  BigReal y,
1175  BigReal z,
1176  int atom1,
1177  int atom2,
1178  BigReal* pairLJEnergy,
1179  BigReal* pairGaussEnergy) const
1180 {
1181  //Initialize return energies to zero
1182  *pairLJEnergy = 0.0;
1183  *pairGaussEnergy = 0.0;
1184 
1185  // Linear search for LJ data
1186  int LJIndex = -1;
1187  int LJbegin = pointerToLJBeg[atom1];
1188  int LJend = pointerToLJEnd[atom1];
1189  for(int i = LJbegin; i <= LJend; i++) {
1190  if(indxLJB[i] == atom2) {
1191  LJIndex = i;
1192  break;
1193  }
1194  }
1195 
1196  // Linear search for Gaussian data
1197  int GaussIndex = -1;
1198  int Gaussbegin = pointerToGaussBeg[atom1];
1199  int Gaussend = pointerToGaussEnd[atom1];
1200  for(int i = Gaussbegin; i <= Gaussend; i++) {
1201  if(indxGaussB[i] == atom2) {
1202  GaussIndex = i;
1203  break;
1204  }
1205  }
1206 
1207  if( LJIndex == -1 && GaussIndex == -1) {
1208  return 0;
1209  } else {
1210  // Code to calculate distances because the pair was found in one of the lists
1211  BigReal r2 = x*x + y*y + z*z;
1212  BigReal r = sqrt(r2);
1213  BigReal ri = 1/r;
1214  BigReal ri6 = (ri*ri*ri*ri*ri*ri);
1215  BigReal ri12 = ri6*ri6;
1216  BigReal ri13 = ri12*ri;
1217  BigReal LJ = 0;
1218  BigReal Gauss = 0;
1219  // Code to calculate LJ
1220  if (LJIndex != -1) {
1221  BigReal ri7 = ri6*ri;
1222  LJ = (12*(pairC12[LJIndex]*ri13) - 6*(pairC6[LJIndex]*ri7));
1223  *pairLJEnergy = (pairC12[LJIndex]*ri12 - pairC6[LJIndex]*ri6);
1224  //std::cout << pairC12[LJIndex] << " " << pairC6[LJIndex] << " " << ri13 << " " << ri7 << " " << LJ << " " << r << "\n";
1225  }
1226  // Code to calculate Gaussian
1227  if (GaussIndex != -1) {
1228  BigReal gr = 12*gRepulsive[GaussIndex]*ri13;
1229  BigReal r1prime = r - gMu1[GaussIndex];
1230  BigReal tmp1 = r1prime * r1prime;
1231  BigReal r2prime = r - gMu2[GaussIndex];
1232  BigReal tmp2 = r2prime * r2prime;
1233  BigReal tmp_gauss1 = 0;
1234  BigReal one_gauss1 = 1;
1235  BigReal tmp_gauss2 = 0;
1236  BigReal one_gauss2 = 1;
1237  if (giSigma1[GaussIndex] != 0) {
1238  tmp_gauss1 = exp(-tmp1*giSigma1[GaussIndex]);
1239  one_gauss1 = 1 - tmp_gauss1;
1240  }
1241  if (giSigma2[GaussIndex] != 0) {
1242  tmp_gauss2 = exp(-tmp2*giSigma2[GaussIndex]);
1243  one_gauss2 = 1 - tmp_gauss2;
1244  }
1245  BigReal A = gA[GaussIndex];
1246  Gauss = gr*one_gauss1*one_gauss2 - 2*A*tmp_gauss1*one_gauss2*r1prime*giSigma1[GaussIndex] \
1247  - 2*tmp_gauss1*one_gauss2*r1prime*giSigma1[GaussIndex]*gRepulsive[GaussIndex]*ri12 - 2*A*tmp_gauss2*one_gauss1*r2prime*giSigma2[GaussIndex] \
1248  - 2*tmp_gauss2*one_gauss1*r2prime*giSigma2[GaussIndex]*gRepulsive[GaussIndex]*ri12;
1249  *pairGaussEnergy = A*(-1+(one_gauss1)*(one_gauss2)*(1+gRepulsive[GaussIndex]*ri12/A));
1250  }
1251  //std::cout << "Net force: " << (LJ + Gauss) << " with ri " << (LJ + Gauss)*ri << "\n";
1252  return (LJ + Gauss)*ri;
1253  }
1254  return 0;
1255 }
1256 // End of get_gro_force2
1257 // JLai
1258 
1259 // JE
1260 BigReal Molecule::get_go_force(BigReal r,
1261  int atom1,
1262  int atom2,
1263  BigReal* goNative,
1264  BigReal* goNonnative) const
1265 {
1266  BigReal goForce = 0.0;
1267  Real pow1;
1268  Real pow2;
1269  // determine which Go chain pair we are working with
1270  //DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ")" << std::endl);
1271  int32 chain1 = atomChainTypes[atom1];
1272  int32 chain2 = atomChainTypes[atom2];
1273 
1274  //DebugM(3," chain1:" << chain1 << ", chain2:" << chain2 << std::endl);
1275  if (chain1 == 0 || chain2 == 0) return 0.0;
1276 
1277  // retrieve Go cutoff for this chain pair
1278  //TMP// JLai -- I'm going to replace this with a constant accessor. This is just a temporary thing
1279  Real goCutoff = const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2);
1280  //DebugM(3," goCutoff:" << goCutoff << std::endl);
1281  if (goCutoff == 0) return 0.0;
1282  // if repulsive then calculate repulsive
1283  // sigmas are initially set to -1.0 if the atom pair fails go_restricted
1284  if (goSigmas[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]] != -1.0) {
1285  if (!goWithinCutoff[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]]) {
1286  Real epsilonRep = const_cast<Molecule*>(this)->get_go_epsilonRep(chain1,chain2);
1287  Real sigmaRep = const_cast<Molecule*>(this)->get_go_sigmaRep(chain1,chain2);
1288  int exp_rep = const_cast<Molecule*>(this)->get_go_exp_rep(chain1,chain2);
1289  pow1 = pow(sigmaRep/r,exp_rep);
1290  goForce = 4*((exp_rep/(r*r)) * epsilonRep * pow1);
1291  *goNative = 0.0;
1292  *goNonnative = (4 * epsilonRep * pow1 );
1293  //DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", epsilonRep:" << epsilonRep << ", sigmaRep:" << sigmaRep << ", r:" << r << ", goForce:" << goForce << std::endl);
1294  }
1295  // if attractive then calculate attractive
1296  else {
1297  int goSigmaIndex1 = goSigmaIndices[atom1];
1298  int goSigmaIndex2 = goSigmaIndices[atom2];
1299  if (goSigmaIndex1 != -1 && goSigmaIndex2 != -1) {
1300  Real epsilon = const_cast<Molecule*>(this)->get_go_epsilon(chain1,chain2);
1301  int exp_a = const_cast<Molecule*>(this)->get_go_exp_a(chain1,chain2);
1302  int exp_b = const_cast<Molecule*>(this)->get_go_exp_b(chain1,chain2);
1303  Real sigma_ij = goSigmas[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]];
1304  // Positive gradient of potential, not negative gradient of potential
1305  pow1 = pow(sigma_ij/r,exp_a);
1306  pow2 = pow(sigma_ij/r,exp_b);
1307  goForce = ((4/(r*r)) * epsilon * (exp_a * pow1 - exp_b * pow2));
1308  //DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", sigma_ij:" << sigma_ij << ", r:" << r << ", goForce:" << goForce << std::endl);
1309  *goNative = (4 * epsilon * ( pow1 - pow2 ) );
1310  *goNonnative = 0.0;
1311  }
1312  }
1313  }
1314  //DebugM(3,"goForce:" << goForce << std::endl);
1315  return goForce;
1316 }
1317 /* END OF FUNCTION get_go_force_old */
1318 
1319 
1320  /************************************************************************/
1321  /* */
1322  /* JE - FUNCTION get_go_force_new */
1323  /* */
1324  /* INPUTS: */
1325  /* r - distance between the two atoms */
1326  /* atom1 - the ID of the first atom */
1327  /* atom2 - the ID of the second atom */
1328  /* */
1329  /* This function calculates the Go force between two atoms. If the */
1330  /* atoms do not have Go parameters or sigmas, 0 is returned. */
1331  /* */
1332  /************************************************************************/
1333 // JE
1334 BigReal Molecule::get_go_force_new(BigReal r,
1335  int atom1,
1336  int atom2,
1337  BigReal* goNative,
1338  BigReal* goNonnative) const
1339 {
1340  int resid1;
1341  int resid2;
1342  int residDiff;
1343  Real xcoorDiff;
1344  Real ycoorDiff;
1345  Real zcoorDiff;
1346  Real atomAtomDist;
1347  Real exp_a;
1348  Real exp_b;
1349  Real sigma_ij;
1350  Real epsilon;
1351  Real epsilonRep;
1352  Real sigmaRep;
1353  Real expRep;
1354  Real pow1;
1355  Real pow2;
1356 
1357  BigReal goForce = 0.0;
1358  *goNative = 0;
1359  *goNonnative = 0;
1360 
1361  // determine which Go chain pair we are working with
1362  DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ")" << std::endl);
1363  int goIndex1 = goSigmaIndices[atom1];
1364  int goIndex2 = goSigmaIndices[atom2];
1365 
1366  int32 chain1 = atomChainTypes[goIndex1];
1367  int32 chain2 = atomChainTypes[goIndex2];
1368 
1369  DebugM(3," chain1:" << chain1 << ", chain2:" << chain2 << std::endl);
1370  if (chain1 == 0 || chain2 == 0) return 0.0;
1371 
1372  // retrieve Go cutoff for this chain pair
1373  Real goCutoff = const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2);
1374  DebugM(3," goCutoff:" << goCutoff << std::endl);
1375  if (goCutoff == 0) return 0.0;
1376 
1377  // sigmas are initially set to -1.0 if the atom pair fails go_restricted
1378  // no goSigmas array anymore
1379  //Real sigma_ij = goSigmas[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]];
1380 
1381  // XXX - used to be a condition for the following if
1382  //if the atoms are within 4 of each other
1383  //if ( !atoms_1to4(atom1,atom2) ) {
1384 
1385  // if goSigmaIndices aren't defined, don't calculate forces
1386  if ( goIndex1 != -1 && goIndex2 != -1 ) {
1387  resid1 = goResids[goIndex1];
1388  resid2 = goResids[goIndex2];
1389  residDiff = resid2 - resid1;
1390  if (residDiff < 0) residDiff = -residDiff;
1391  // if this is a Go atom pair that is not restricted
1392  if ( !(const_cast<Molecule*>(this)->go_restricted(chain1,chain2,residDiff)) ) {
1393  xcoorDiff = goCoordinates[goIndex1*3] - goCoordinates[goIndex2*3];
1394  ycoorDiff = goCoordinates[goIndex1*3 + 1] - goCoordinates[goIndex2*3 + 1];
1395  zcoorDiff = goCoordinates[goIndex1*3 + 2] - goCoordinates[goIndex2*3 + 2];
1396  atomAtomDist = sqrt(xcoorDiff*xcoorDiff + ycoorDiff*ycoorDiff + zcoorDiff*zcoorDiff);
1397 
1398  // if attractive then calculate attractive
1399  if ( atomAtomDist <= const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2) ) {
1400  exp_a = const_cast<Molecule*>(this)->get_go_exp_a(chain1,chain2);
1401  exp_b = const_cast<Molecule*>(this)->get_go_exp_b(chain1,chain2);
1402  sigma_ij = pow(static_cast<double>(exp_b/exp_a),(1.0/(exp_a-exp_b))) * atomAtomDist;
1403 
1404  // [JLai] print out atoms involved in native contacts
1405  // printf("ATOM1: %d C1: %d ATOM2: %d C2: %d\n", atom1,chain1,atom2,chain2);
1406 
1407  epsilon = const_cast<Molecule*>(this)->get_go_epsilon(chain1,chain2);
1408  pow1 = pow(sigma_ij/r,static_cast<double>(exp_a));
1409  pow2 = pow(sigma_ij/r,static_cast<double>(exp_b));
1410  //goForce = ((4/r) * epsilon * (exp_a * pow(sigma_ij/r,exp_a) - exp_b * pow(sigma_ij/r,exp_b)));
1411  goForce = ((4/(r*r)) * epsilon * (exp_a * pow1 - exp_b * pow2));
1412  DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", exp_a:" << exp_a << ", exp_b:" << exp_b << ", sigma_ij:" << sigma_ij << ", r:" << r << ", goForce:" << goForce << std::endl);
1413  //goEnergy = (4 * epsilon * ( pow(sigma_ij/r,exp_a) - pow(sigma_ij/r,exp_b) ) ); // JLai I changed some of the expressions
1414  *goNative = (4 * epsilon * ( pow1 - pow2 ) );
1415  *goNonnative = 0;
1416  }
1417 
1418  // if repulsive then calculate repulsive
1419  else {
1420  epsilonRep = const_cast<Molecule*>(this)->get_go_epsilonRep(chain1,chain2);
1421  sigmaRep = const_cast<Molecule*>(this)->get_go_sigmaRep(chain1,chain2);
1422  expRep = const_cast<Molecule*>(this)->get_go_exp_rep(chain1,chain2);
1423  pow1 = pow(sigmaRep/r,(BigReal)expRep);
1424  //goForce = ((12.0/r) * epsilonRep * pow(sigmaRep/r,12.0));
1425  goForce = ((4/(r*r)) * expRep * epsilonRep * pow1);
1426  DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", epsilonRep:" << epsilonRep << ", sigmaRep:" << sigmaRep << ", r:" << r << ", goForce:" << goForce << std::endl);
1427  //goEnergy = (4 * epsilonRep * pow(sigmaRep/r,12.0)); // JLai I changed some of the expressions
1428  *goNonnative = (4 * epsilonRep * pow1);
1429  *goNative = 0;
1430  }
1431  }
1432  }
1433 
1434  //DebugM(3,"goForce:" << goForce << std::endl);
1435  return goForce;
1436 }
1437 /* END OF FUNCTION get_go_force_new */
1438 
1439 
1440  /************************************************************************/
1441  /* */
1442  /* JLai - FUNCTION get_go_force2 */
1443  /* */
1444  /* INPUTS: */
1445  /* x - the x distance between p_i and p_j */
1446  /* y - the y distance between p_i and p_j */
1447  /* z - the z distance between p_i and p_j */
1448  /* atom1 - the ID of the second atom */
1449  /* atom2 - the ID of the second atom */
1450  /* */
1451  /* This function returns the force between two input atoms give their */
1452 /* distance and their atom indices. */
1453  /* */
1454  /************************************************************************/
1455 // JLai
1456 BigReal Molecule::get_go_force2(BigReal x,
1457  BigReal y,
1458  BigReal z,
1459  int atom1,
1460  int atom2,
1461  BigReal *goNative,
1462  BigReal *goNonnative) const
1463 {
1464 
1465  // Check to see if restricted. If so, escape function early
1466  int32 chain1 = atomChainTypes[atom1];
1467  int32 chain2 = atomChainTypes[atom2];
1468 
1469  if(chain1 == 0 || chain2 == 0) return 0.0;
1470  Molecule *mol = const_cast<Molecule*>(this);
1471  Real goCutoff = mol->get_go_cutoff(chain1,chain2);
1472  if(goCutoff == 0) return 0.0;
1473 
1474  int resid1 = goResidIndices[atom1];
1475  int resid2 = goResidIndices[atom2];
1476  int residDiff = abs(resid1 - resid2);
1477  if((mol->go_restricted(chain1,chain2,residDiff))) {
1478  return 0.0;
1479  }
1480 
1481  int LJIndex = -1;
1482  int LJbegin = pointerToGoBeg[atom1];
1483  int LJend = pointerToGoEnd[atom1];
1484  for(int i = LJbegin; i <= LJend; i++) {
1485  if(goIndxLJB[i] == atom2) {
1486  LJIndex = i;
1487  }
1488  }
1489 
1490  BigReal r2 = x*x + y*y + z*z;
1491  BigReal r = sqrt(r2);
1492 
1493  if (LJIndex == -1) {
1494  int exp_rep = const_cast<Molecule*>(this)->get_go_exp_rep(chain1,chain2);
1495  BigReal epsilonRep = const_cast<Molecule*>(this)->get_go_epsilonRep(chain1, chain2);
1496  BigReal sigmaRep = const_cast<Molecule*>(this)->get_go_sigmaRep(chain1, chain2);
1497  double sigmaRepPow = pow(sigmaRep,exp_rep);
1498  BigReal LJ = (4*epsilonRep*exp_rep*sigmaRepPow*pow(r,-(exp_rep+1)));
1499  *goNative = 0;
1500  *goNonnative = (4*epsilonRep*sigmaRepPow*pow(r,-exp_rep));
1501  //*goNonnative = (4*epsilonRep * pow(sigmaRep/r,exp_rep));
1502  return (LJ/r);
1503  } else {
1504  // Code to calculate distances because the pair was found in one of the lists
1505  int exp_a = const_cast<Molecule*>(this)->get_go_exp_a(chain1,chain2);
1506  int exp_b = const_cast<Molecule*>(this)->get_go_exp_b(chain1,chain2);
1507  // We want the force, so we have to take the n+1 derivative
1508  BigReal powA = pow(r,-(exp_a + 1));
1509  BigReal powB = pow(r,-(exp_b + 1));
1510  BigReal powaa = pow(r,-exp_a);
1511  BigReal powbb = pow(r,-exp_b);
1512  BigReal epsilon = const_cast<Molecule*>(this)->get_go_epsilon(chain1,chain2);
1513  BigReal LJ = 4 * epsilon * (exp_a*goSigmaPairA[LJIndex]*powA - exp_b*goSigmaPairB[LJIndex]*powB);
1514  *goNative = 4 * epsilon * (goSigmaPairA[LJIndex]*powaa - goSigmaPairB[LJIndex]*powbb);
1515  *goNonnative = 0;
1516  return (LJ/r);
1517  }
1518  return 0;
1519 }
1520 // JLai
1521 /* END OF FUNCTION get_go_force2 */
1522 
1523 #ifndef MEM_OPT_VERSION
1524  /************************************************************************/
1525  /* */
1526  /* JE - FUNCTION atoms_1to4 */
1527  /* */
1528  /* INPUTS: */
1529  /* atom1 - the ID of the first atom */
1530  /* atom2 - the ID of the second atom */
1531  /* */
1532  /* This function tells whether or not the two input atoms are within */
1533  /* 1 to 4 bonds away from each other: bonds, angles, or dihedrals. */
1534  /* */
1535  /************************************************************************/
1536 // JE
1537 Bool Molecule::atoms_1to4(unsigned int atom1,
1538  unsigned int atom2)
1539 {
1540  int bondNum; // Bonds in bonded list
1541  int angleNum; // Angles in angle list
1542  int dihedralNum; // Dihedrals in dihedral list
1543  int *bonds;
1544  int *angles;
1545  int *dihedrals;
1546  Bond *bond; // Temporary bond structure
1547  Angle *angle; // Temporary angle structure
1548  Dihedral *dihedral; // Temporary dihedral structure
1549 
1550  DebugM(2,"atoms_1to4(" << atom1 << "," << atom2 << ")" << std::endl);
1551 
1552  bonds = get_bonds_for_atom(atom1);
1553  bondNum = *bonds;
1554  while(bondNum != -1) {
1555  bond = get_bond(bondNum);
1556  DebugM(2,"bond atom1:" << bond->atom1 << ", atom2:" << bond->atom2 << std::endl);
1557  if (atom2 == bond->atom1 || atom2 == bond->atom2) {
1558  return TRUE;
1559  }
1560  bondNum = *(++bonds);
1561  }
1562 
1563  bonds = get_bonds_for_atom(atom2);
1564  bondNum = *bonds;
1565  while(bondNum != -1) {
1566  bond = get_bond(bondNum);
1567  DebugM(2,"bond atom1:" << bond->atom1 << ", atom2:" << bond->atom2 << std::endl);
1568  if (atom1 == bond->atom1 || atom1 == bond->atom2) {
1569  return TRUE;
1570  }
1571  bondNum = *(++bonds);
1572  }
1573 
1574  angles = get_angles_for_atom(atom1);
1575  angleNum = *angles;
1576  while(angleNum != -1) {
1577  angle = get_angle(angleNum);
1578  DebugM(2,"angle atom1:" << angle->atom1 << ", atom2:" << angle->atom2 << ", atom3:" << angle->atom3 << std::endl);
1579  if (atom2 == angle->atom1 || atom2 == angle->atom2 || atom2 == angle->atom3) {
1580  return TRUE;
1581  }
1582  angleNum = *(++angles);
1583  }
1584 
1585  angles = get_angles_for_atom(atom2);
1586  angleNum = *angles;
1587  while(angleNum != -1) {
1588  angle = get_angle(angleNum);
1589  DebugM(2,"angle atom1:" << angle->atom1 << ", atom2:" << angle->atom2 << ", atom3:" << angle->atom3 << std::endl);
1590  if (atom1 == angle->atom1 || atom1 == angle->atom2 || atom1 == angle->atom3) {
1591  return TRUE;
1592  }
1593  angleNum = *(++angles);
1594  }
1595 
1596  dihedrals = get_dihedrals_for_atom(atom1);
1597  dihedralNum = *dihedrals;
1598  while(dihedralNum != -1) {
1599  dihedral = get_dihedral(dihedralNum);
1600  DebugM(2,"dihedral atom1:" << dihedral->atom1 << ", atom2:" << dihedral->atom2 << ", atom3:" << dihedral->atom3 << ", atom4:" << dihedral->atom4 << std::endl);
1601  if (atom2 == dihedral->atom1 || atom2 == dihedral->atom2 \
1602  || atom2 == dihedral->atom3 || atom2 == dihedral->atom4) {
1603  return TRUE;
1604  }
1605  dihedralNum = *(++dihedrals);
1606  }
1607 
1608  dihedrals = get_dihedrals_for_atom(atom2);
1609  dihedralNum = *dihedrals;
1610  while(dihedralNum != -1) {
1611  dihedral = get_dihedral(dihedralNum);
1612  DebugM(2,"dihedral atom1:" << dihedral->atom1 << ", atom2:" << dihedral->atom2 << ", atom3:" << dihedral->atom3 << ", atom4:" << dihedral->atom4 << std::endl);
1613  if (atom1 == dihedral->atom1 || atom1 == dihedral->atom2 \
1614  || atom1 == dihedral->atom3 || atom1 == dihedral->atom4) {
1615  return TRUE;
1616  }
1617  dihedralNum = *(++dihedrals);
1618  }
1619 
1620  return FALSE;
1621 }
1622 /* END OF FUNCTION atoms_1to4 */
1623 #endif // #ifndef MEM_OPT_VERSION
1624 
1625 //JLai
1626 /************************************************************************/
1627 /* */
1628 /* FUNCTION send_GoMolecule */
1629 /* */
1630 /* send_Molecule is used by the Master node to distribute the */
1631 /* Go information to all the client nodes. It is NEVER called*/
1632 /* by the client nodes. */
1633 /* */
1634 /************************************************************************/
1635 void Molecule::send_GoMolecule(MOStream *msg) {
1636  Real *a1, *a2, *a3, *a4;
1637  int *i1, *i2, *i3, *i4;
1638  int maxGoChainsSqr = MAX_GO_CHAINS*MAX_GO_CHAINS; // JE JLai Go code
1639  msg->put(NumGoChains);
1640 
1641  if (NumGoChains) {
1642  // int go_indices[MAX_GO_CHAINS+1]; // Indices from chainIDs to go_array
1643  // GoValue go_array[MAX_GO_CHAINS*MAX_GO_CHAINS]; // Array of Go params
1644  msg->put(MAX_GO_CHAINS+1,go_indices);
1645 
1646  a1 = new Real[maxGoChainsSqr];
1647  a2 = new Real[maxGoChainsSqr];
1648  a3 = new Real[maxGoChainsSqr];
1649  a4 = new Real[maxGoChainsSqr];
1650  i1 = new int[maxGoChainsSqr];
1651  i2 = new int[maxGoChainsSqr];
1652  i3 = new int[maxGoChainsSqr];
1653  i4 = new int[maxGoChainsSqr*MAX_RESTRICTIONS];
1654 
1655  if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) ||
1656  (i1 == NULL) || (i2 == NULL) || (i3 == NULL) || (i4 == NULL) )
1657  {
1658  NAMD_die("memory allocation failed in Molecules::send_Molecules");
1659  }
1660 
1661  for (int i=0; i<maxGoChainsSqr; i++) {
1662  a1[i] = go_array[i].epsilon;
1663  a2[i] = go_array[i].sigmaRep;
1664  a3[i] = go_array[i].epsilonRep;
1665  a4[i] = go_array[i].cutoff;
1666  i1[i] = go_array[i].exp_a;
1667  i2[i] = go_array[i].exp_b;
1668  i3[i] = go_array[i].exp_rep;
1669  for (int j=0; j<MAX_RESTRICTIONS; j++) {
1670  i4[i*MAX_RESTRICTIONS + j] = go_array[i].restrictions[j];
1671  }
1672  }
1673 
1674  msg->put(maxGoChainsSqr, a1);
1675  msg->put(maxGoChainsSqr, a2);
1676  msg->put(maxGoChainsSqr, a3);
1677  msg->put(maxGoChainsSqr, a4);
1678  msg->put(maxGoChainsSqr, i1);
1679  msg->put(maxGoChainsSqr, i2);
1680  msg->put(maxGoChainsSqr, i3);
1681  msg->put(maxGoChainsSqr*MAX_RESTRICTIONS, i4);
1682 
1683  delete [] a1;
1684  delete [] a2;
1685  delete [] a3;
1686  delete [] a4;
1687  delete [] i1;
1688  delete [] i2;
1689  delete [] i3;
1690  delete [] i4;
1691  }
1692 
1693  //Ported JLai
1694  if (simParams->goForcesOn) {
1695  switch(simParams->goMethod) {
1696  case 1:
1697  msg->put(numGoAtoms);
1698  msg->put(numAtoms, goSigmaIndices);
1699  msg->put(numGoAtoms, atomChainTypes);
1700  msg->put(numGoAtoms*numGoAtoms, goSigmas);
1701  msg->put(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
1702  // printf("Molecule.C sending atomChainTypes %d %d \n", numGoAtoms, atomChainTypes);
1703  break;
1704  case 2: //GSS
1705  msg->put(numGoAtoms);
1706  msg->put(numAtoms,pointerToGoBeg);
1707  msg->put(numAtoms,pointerToGoEnd);
1708  msg->put(numAtoms,goSigmaIndices);
1709  msg->put(numAtoms,goResidIndices);
1710  msg->put(numGoAtoms,atomChainTypes);
1711  msg->put(goNumLJPair);
1712  msg->put(goNumLJPair,goIndxLJA);
1713  msg->put(goNumLJPair,goIndxLJB);
1714  msg->put(goNumLJPair,goSigmaPairA);
1715  msg->put(goNumLJPair,goSigmaPairB);
1716  break;
1717  case 3:
1718  msg->put(numGoAtoms);
1719  msg->put(numAtoms, goSigmaIndices);
1720  msg->put(numGoAtoms, atomChainTypes);
1721  //msg->put(numGoAtoms*numGoAtoms, goSigmas);
1722  //msg->put(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
1723  msg->put(numGoAtoms*3, goCoordinates);
1724  msg->put(numGoAtoms, goResids);
1725  break;
1726  }
1727  }
1728 
1729  msg->end();
1730  delete msg;
1731 }
1732 /* END OF FUNCTION send_GoMolecule */
1733 
1734 // JLai
1735 /************************************************************************/
1736 /* */
1737 /* FUNCTION receive_Molecule */
1738 /* */
1739 /* receive_Molecule is used by all the clients to receive the */
1740 /* Go structural data sent out by the master node. It is NEVER called */
1741 /* by the Master node. */
1742 /* */
1743 /************************************************************************/
1744 void Molecule::receive_GoMolecule(MIStream *msg) {
1745  // Ported by JLai -- Original by JE
1746  // JE - receive Go info
1747  Real *a1, *a2, *a3, *a4;
1748  int *i1, *i2, *i3, *i4;
1749  int maxGoChainsSqr = MAX_GO_CHAINS*MAX_GO_CHAINS; // JE JLai Go code
1750  msg->get(NumGoChains);
1751 
1752  if (NumGoChains) {
1753  //go_indices = new int[MAX_GO_CHAINS+1];
1754  //go_array = new GoValue[MAX_GO_CHAINS*MAX_GO_CHAINS];
1755 
1756  // int go_indices[MAX_GO_CHAINS+1]; // Indices from chainIDs to go_array
1757  // GoValue go_array[MAX_GO_CHAINS*MAX_GO_CHAINS]; // Array of Go params
1758  msg->get(MAX_GO_CHAINS+1,go_indices);
1759 
1760  a1 = new Real[maxGoChainsSqr];
1761  a2 = new Real[maxGoChainsSqr];
1762  a3 = new Real[maxGoChainsSqr];
1763  a4 = new Real[maxGoChainsSqr];
1764  i1 = new int[maxGoChainsSqr];
1765  i2 = new int[maxGoChainsSqr];
1766  i3 = new int[maxGoChainsSqr];
1767  i4 = new int[maxGoChainsSqr*MAX_RESTRICTIONS];
1768 
1769  if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) ||
1770  (i1 == NULL) || (i2 == NULL) || (i3 == NULL) || (i4 == NULL) )
1771  {
1772  NAMD_die("memory allocation failed in Molecule::send_Molecule");
1773  }
1774 
1775  msg->get(maxGoChainsSqr, a1);
1776  msg->get(maxGoChainsSqr, a2);
1777  msg->get(maxGoChainsSqr, a3);
1778  msg->get(maxGoChainsSqr, a4);
1779  msg->get(maxGoChainsSqr, i1);
1780  msg->get(maxGoChainsSqr, i2);
1781  msg->get(maxGoChainsSqr, i3);
1782  msg->get(maxGoChainsSqr*MAX_RESTRICTIONS, i4);
1783 
1784  for (int i=0; i<maxGoChainsSqr; i++) {
1785  go_array[i].epsilon = a1[i];
1786  go_array[i].sigmaRep = a2[i];
1787  go_array[i].epsilonRep = a3[i];
1788  go_array[i].cutoff = a4[i];
1789  go_array[i].exp_a = i1[i];
1790  go_array[i].exp_b = i2[i];
1791  go_array[i].exp_rep = i3[i];
1792  for (int j=0; j<MAX_RESTRICTIONS; j++) {
1793  go_array[i].restrictions[j] = i4[i*MAX_RESTRICTIONS + j];
1794  }
1795  }
1796 
1797  delete [] a1;
1798  delete [] a2;
1799  delete [] a3;
1800  delete [] a4;
1801  delete [] i1;
1802  delete [] i2;
1803  delete [] i3;
1804  delete [] i4;
1805 
1806  //msg->get(MAX_GO_CHAINS*MAX_GO_CHAINS, (char*)go_array);
1807 
1808  /*DebugM(3,"NumGoChains:" << NumGoChains << std::endl);
1809  for (int ii=0; ii<MAX_GO_CHAINS; ii++) {
1810  for (int jj=0; jj<MAX_GO_CHAINS; jj++) {
1811  DebugM(3,"go_array[" << ii << "][" << jj << "]:" << go_array[ii][jj] << std::endl);
1812  }
1813  }
1814  for (int ii=0; ii<MAX_GO_CHAINS+1; ii++) {
1815  DebugM(3,"go_indices[" << ii << "]:" << go_indices[ii] << std::endl);
1816  }*/
1817  }
1818 
1819  if (simParams->goForcesOn) {
1820  switch(simParams->goMethod) {
1821  case 1:
1822  msg->get(numGoAtoms);
1823  //printf("Deleting goSigmaIndiciesA\n");
1824  delete [] goSigmaIndices;
1825  goSigmaIndices = new int32[numAtoms];
1826  //printf("Deleting atomChainTypesA\n");
1827  delete [] atomChainTypes;
1828  atomChainTypes = new int32[numGoAtoms];
1829  //printf("Deleting goSigmasA\n");
1830  delete [] goSigmas;
1831  goSigmas = new Real[numGoAtoms*numGoAtoms];
1832  //printf("Deleting goWithinCutoffA\n");
1833  delete [] goWithinCutoff;
1834  goWithinCutoff = new bool[numGoAtoms*numGoAtoms];
1835  msg->get(numAtoms, goSigmaIndices);
1836  msg->get(numGoAtoms, atomChainTypes);
1837  msg->get(numGoAtoms*numGoAtoms, goSigmas);
1838  msg->get(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
1839  break;
1840  case 2: //GSR
1841  msg->get(numGoAtoms);
1842  delete [] pointerToGoBeg;
1843  pointerToGoBeg = new int[numAtoms];
1844  msg->get(numAtoms,pointerToGoBeg);
1845  delete [] pointerToGoEnd;
1846  pointerToGoEnd = new int[numAtoms];
1847  msg->get(numAtoms,pointerToGoEnd);
1848  delete [] goSigmaIndices;
1849  goSigmaIndices = new int32[numAtoms];
1850  msg->get(numAtoms,goSigmaIndices);
1851  delete [] goResidIndices;
1852  goResidIndices = new int32[numAtoms];
1853  msg->get(numAtoms,goResidIndices);
1854  delete [] atomChainTypes;
1855  atomChainTypes = new int32[numGoAtoms];
1856  msg->get(numGoAtoms,atomChainTypes);
1857  msg->get(goNumLJPair);
1858  delete [] goIndxLJA;
1859  goIndxLJA = new int[goNumLJPair];
1860  msg->get(goNumLJPair,goIndxLJA);
1861  delete [] goIndxLJB;
1862  goIndxLJB = new int[goNumLJPair];
1863  msg->get(goNumLJPair,goIndxLJB);
1864  delete [] goSigmaPairA;
1865  goSigmaPairA = new double[goNumLJPair];
1866  msg->get(goNumLJPair,goSigmaPairA);
1867  delete [] goSigmaPairB;
1868  goSigmaPairB = new double[goNumLJPair];
1869  msg->get(goNumLJPair,goSigmaPairB);
1870  break;
1871  case 3:
1872  msg->get(numGoAtoms);
1873  //printf("Deleting goSigmaIndiciesB\n");
1874  delete [] goSigmaIndices;
1875  goSigmaIndices = new int32[numAtoms];
1876  //printf("Deleting atomChainTypesB\n");
1877  delete [] atomChainTypes;
1878  atomChainTypes = new int32[numGoAtoms];
1879  //delete [] goSigmas;
1880  //goSigmas = new Real[numGoAtoms*numGoAtoms];
1881  //delete [] goWithinCutoff;
1882  //goWithinCutoff = new bool[numGoAtoms*numGoAtoms];
1883  //printf("Deleting goCoordinatesB\n");
1884  delete [] goCoordinates;
1885  goCoordinates = new Real[numGoAtoms*3];
1886  //printf("Deleting goResidsB\n");
1887  delete [] goResids;
1888  goResids = new int[numGoAtoms];
1889  msg->get(numAtoms, goSigmaIndices);
1890  msg->get(numGoAtoms, atomChainTypes);
1891  //msg->get(numGoAtoms*numGoAtoms, goSigmas);
1892  //msg->get(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
1893  msg->get(numGoAtoms*3, goCoordinates);
1894  msg->get(numGoAtoms, goResids);
1895  break;
1896  }
1897  }
1898 
1899  delete msg;
1900 
1901 }
1902 /* END OF FUNCTION receive_GoMolecule */
Molecule::print_go_params
void print_go_params()
Definition: GoMolecule.C:548
go_val::exp_b
int exp_b
Definition: Molecule.h:108
Molecule::get_go_exp_a
int get_go_exp_a(int chain1, int chain2)
Definition: Molecule.h:1613
iINFO
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
Molecule::get_go_cutoff
Real get_go_cutoff(int chain1, int chain2)
Definition: Molecule.h:1601
Molecule::get_dihedrals_for_atom
int32 * get_dihedrals_for_atom(int anum)
Definition: Molecule.h:1149
MOStream::end
void end(void)
Definition: MStream.C:176
PDB
Definition: PDB.h:35
Molecule::gMu2
Real * gMu2
Definition: Molecule.h:679
go_pair::A
double A
Definition: Molecule.h:120
NAMD_read_line
int NAMD_read_line(FILE *fd, char *buf, int bufsize)
Definition: strlib.C:38
Molecule::send_GoMolecule
void send_GoMolecule(MOStream *)
Definition: GoMolecule.C:1635
Molecule::pointerToGoBeg
int * pointerToGoBeg
Definition: Molecule.h:657
Molecule::giSigma2
Real * giSigma2
Definition: Molecule.h:680
Molecule::goPDB
PDB * goPDB
Definition: Molecule.h:650
int32
short int32
Definition: dumpdcd.c:24
Molecule::get_go_force_new
BigReal get_go_force_new(BigReal, int, int, BigReal *, BigReal *) const
Definition: GoMolecule.C:1334
go_val::epsilonRep
Real epsilonRep
Definition: Molecule.h:111
Molecule::get_gro_force2
BigReal get_gro_force2(BigReal, BigReal, BigReal, int, int, BigReal *, BigReal *) const
Definition: GoMolecule.C:1173
Molecule::goCoordinates
Real * goCoordinates
Definition: Molecule.h:648
Molecule::receive_GoMolecule
void receive_GoMolecule(MIStream *)
Definition: GoMolecule.C:1744
dihedral::atom3
int32 atom3
Definition: structures.h:65
MAX_RESTRICTIONS
#define MAX_RESTRICTIONS
Definition: Molecule.h:37
Molecule::get_go_epsilonRep
Real get_go_epsilonRep(int chain1, int chain2)
Definition: Molecule.h:1604
A
const BigReal A
Definition: ComputeNonbondedBase2.h:209
Molecule::pairC6
Real * pairC6
Definition: Molecule.h:667
go_pair::goIndxB
int goIndxB
Definition: Molecule.h:119
Real
float Real
Definition: common.h:109
DebugM
#define DebugM(x, y)
Definition: Debug.h:59
Molecule::goSigmaPairA
double * goSigmaPairA
Definition: Molecule.h:655
FALSE
#define FALSE
Definition: common.h:118
Molecule::pointerToLJEnd
int * pointerToLJEnd
Definition: Molecule.h:664
iWARN
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
MIStream::get
MIStream * get(char &data)
Definition: MStream.h:29
go_pair
Definition: Molecule.h:116
Molecule::goWithinCutoff
bool * goWithinCutoff
Definition: Molecule.h:647
PDBAtom::zcoor
BigReal zcoor(void)
Definition: PDBData.C:433
NAMD_find_first_word
void NAMD_find_first_word(char *source, char *word)
Definition: strlib.C:258
iout
#define iout
Definition: InfoStream.h:87
NAMD_FILENAME_BUFFER_SIZE
#define NAMD_FILENAME_BUFFER_SIZE
Definition: common.h:38
PDB::num_atoms
int num_atoms(void)
Definition: PDB.C:323
Molecule::gA
Real * gA
Definition: Molecule.h:676
MAX_GO_CHAINS
#define MAX_GO_CHAINS
Definition: Molecule.h:36
Molecule::read_go_file
void read_go_file(char *)
Definition: GoMolecule.C:113
Molecule::build_go_arrays
void build_go_arrays(StringList *, char *)
Definition: GoMolecule.C:950
go_val::exp_a
int exp_a
Definition: Molecule.h:107
Molecule::get_angles_for_atom
int32 * get_angles_for_atom(int anum)
Definition: Molecule.h:1147
Molecule::indxLJB
int * indxLJB
Definition: Molecule.h:666
Molecule::get_go_epsilon
Real get_go_epsilon(int chain1, int chain2)
Definition: Molecule.h:1610
Molecule::pointerToGaussBeg
int * pointerToGaussBeg
Definition: Molecule.h:671
Molecule::atomChainTypes
int32 * atomChainTypes
Definition: Molecule.h:643
dihedral::atom4
int32 atom4
Definition: structures.h:66
bond::atom1
int32 atom1
Definition: structures.h:48
Molecule::gMu1
Real * gMu1
Definition: Molecule.h:677
go_pair::B
double B
Definition: Molecule.h:121
Molecule::pointerToLJBeg
int * pointerToLJBeg
Definition: Molecule.h:663
Molecule::gRepulsive
Real * gRepulsive
Definition: Molecule.h:681
Molecule::goInit
void goInit()
Definition: GoMolecule.C:54
go_pair::goIndxA
int goIndxA
Definition: Molecule.h:118
Molecule::pointerToGoEnd
int * pointerToGoEnd
Definition: Molecule.h:658
Molecule::indxGaussB
int * indxGaussB
Definition: Molecule.h:675
Molecule::get_go_sigmaRep
Real get_go_sigmaRep(int chain1, int chain2)
Definition: Molecule.h:1607
Molecule::get_go_force
BigReal get_go_force(BigReal, int, int, BigReal *, BigReal *) const
Definition: GoMolecule.C:1260
Molecule::goIndxLJA
int * goIndxLJA
Definition: Molecule.h:653
go_val
Definition: Molecule.h:104
Molecule::get_dihedral
Dihedral * get_dihedral(int dnum) const
Definition: Molecule.h:1074
PDB::atom
PDBAtom * atom(int place)
Definition: PDB.C:394
Molecule::goNumLJPair
int goNumLJPair
Definition: Molecule.h:652
NAMD_blank_string
int NAMD_blank_string(char *str)
Definition: strlib.C:222
Molecule::get_go_exp_rep
int get_go_exp_rep(int chain1, int chain2)
Definition: Molecule.h:1619
angle::atom2
int32 atom2
Definition: structures.h:56
ResizeArray::end
iterator end(void)
Definition: ResizeArray.h:37
PDBAtom::ycoor
BigReal ycoor(void)
Definition: PDBData.C:429
SimParameters::goMethod
GoChoices goMethod
Definition: SimParameters.h:304
Molecule::goResidIndices
int32 * goResidIndices
Definition: Molecule.h:645
Molecule::go_array
GoValue go_array[MAX_GO_CHAINS *MAX_GO_CHAINS]
Definition: Molecule.h:1584
z
gridSize z
Definition: ComputePmeCUDAKernel.cu:427
angle::atom1
int32 atom1
Definition: structures.h:55
Bool
int Bool
Definition: common.h:133
Molecule::pointerToGaussEnd
int * pointerToGaussEnd
Definition: Molecule.h:672
PDBAtom::residueseq
int residueseq(void)
Definition: PDBData.C:411
Molecule::get_bonds_for_atom
int32 * get_bonds_for_atom(int anum)
Definition: Molecule.h:1145
Molecule::goResids
int * goResids
Definition: Molecule.h:649
Molecule::goSigmas
Real * goSigmas
Definition: Molecule.h:646
angle::atom3
int32 atom3
Definition: structures.h:57
Molecule::energyNonnative
BigReal energyNonnative
Definition: Molecule.h:685
Molecule::go_restricted
Bool go_restricted(int, int, int)
Definition: GoMolecule.C:525
go_val::cutoff
Real cutoff
Definition: Molecule.h:112
Molecule::get_bond
Bond * get_bond(int bnum) const
Definition: Molecule.h:1065
Molecule::pairC12
Real * pairC12
Definition: Molecule.h:668
Molecule::numAtoms
int numAtoms
Definition: Molecule.h:556
go_val::sigmaRep
Real sigmaRep
Definition: Molecule.h:110
NAMD_die
void NAMD_die(const char *err_msg)
Definition: common.C:83
goNonnative
BigReal goNonnative
Definition: ComputeNonbondedBase2.h:431
PDBAtom::xcoor
BigReal xcoor(void)
Definition: PDBData.C:425
dihedral::atom2
int32 atom2
Definition: structures.h:64
Molecule::goSigmaPairB
double * goSigmaPairB
Definition: Molecule.h:656
ResizeArray::add
int add(const Elem &elem)
Definition: ResizeArray.h:97
sort
BlockRadixSort::TempStorage sort
Definition: CudaTileListKernel.cu:654
Molecule::get_go_exp_b
int get_go_exp_b(int chain1, int chain2)
Definition: Molecule.h:1616
StringList::next
StringList * next
Definition: ConfigList.h:49
goForce
BigReal goForce
Definition: ComputeNonbondedBase2.h:432
StringList::data
char * data
Definition: ConfigList.h:48
Molecule::get_go_force2
BigReal get_go_force2(BigReal, BigReal, BigReal, int, int, BigReal *, BigReal *) const
Definition: GoMolecule.C:1456
Molecule::go_indices
int go_indices[MAX_GO_CHAINS+1]
Definition: Molecule.h:1585
Molecule::build_go_sigmas
void build_go_sigmas(StringList *, char *)
Definition: GoMolecule.C:577
Molecule::giSigma1
Real * giSigma1
Definition: Molecule.h:678
y
gridSize y
Definition: ComputePmeCUDAKernel.cu:426
MOStream::put
MOStream * put(char data)
Definition: MStream.h:112
Molecule::energyNative
BigReal energyNative
Definition: Molecule.h:684
Molecule::print_go_sigmas
void print_go_sigmas()
Definition: GoMolecule.C:1134
Molecule::NumGoChains
int NumGoChains
Definition: Molecule.h:1586
endi
infostream & endi(infostream &s)
Definition: InfoStream.C:38
go_val::restrictions
int restrictions[MAX_RESTRICTIONS]
Definition: Molecule.h:113
go_val::epsilon
Real epsilon
Definition: Molecule.h:106
Molecule::atoms_1to4
Bool atoms_1to4(unsigned int, unsigned int)
Definition: GoMolecule.C:1537
Molecule::goSigmaIndices
int32 * goSigmaIndices
Definition: Molecule.h:644
x
gridSize x
Definition: ComputePmeCUDAKernel.cu:425
Molecule::goIndxLJB
int * goIndxLJB
Definition: Molecule.h:654
Molecule::get_angle
Angle * get_angle(int anum) const
Definition: Molecule.h:1068
TRUE
#define TRUE
Definition: common.h:119
go_val::exp_rep
int exp_rep
Definition: Molecule.h:109
Molecule::numGoAtoms
int numGoAtoms
Definition: Molecule.h:642
dihedral::atom1
int32 atom1
Definition: structures.h:63
bond::atom2
int32 atom2
Definition: structures.h:49
BigReal
double BigReal
Definition: common.h:114
Molecule::numLJPair
int numLJPair
Definition: Molecule.h:661
Molecule::build_go_params
void build_go_params(StringList *)
Definition: GoMolecule.C:79
Molecule::build_go_sigmas2
void build_go_sigmas2(StringList *, char *)
Definition: GoMolecule.C:747
ResizeArray::begin
iterator begin(void)
Definition: ResizeArray.h:36
Generated by   doxygen 1.8.5