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