GoMolecule.C

Go to the documentation of this file.
00001 
00007 #include <stdio.h>
00008 #include <string.h>
00009 #include <stdlib.h>
00010 #include <ctype.h>
00011 #include "ResizeArray.h"
00012 #include "InfoStream.h"
00013 #include "Molecule.h"
00014 #include "strlib.h"
00015 #include "MStream.h"
00016 #include "Communicate.h"
00017 #include "Node.h"
00018 #include "ObjectArena.h"
00019 #include "Parameters.h"
00020 #include "PDB.h"
00021 #include "SimParameters.h"
00022 #include "Hydrogen.h"
00023 #include "UniqueSetIter.h"
00024 #include "charm++.h"
00025 /* BEGIN gf */
00026 #include "ComputeGridForce.h"
00027 #include "GridForceGrid.h"
00028 #include "MGridforceParams.h"
00029 /* END gf */
00030 
00031 #define MIN_DEBUG_LEVEL 3
00032 //#define DEBUGM
00033 #include "Debug.h"
00034 #include "CompressPsf.h"
00035 #include "ParallelIOMgr.h"
00036 #include <deque>
00037 #include <algorithm>
00038 #ifndef M_PI
00039 #define M_PI 3.14159265358979323846
00040 #endif
00041 #ifndef CODE_REDUNDANT
00042 #define CODE_REDUNDANT 0
00043 #endif
00044 
00045 // Ported by JLai from NAMD 2.7
00046 /************************************************************************/
00047 /*                                                                      */
00048 /*      FUNCTION goInit                                                 */
00049 /*                                                                      */
00050 /*        This function is only called from the Molecule constructor.   */
00051 /*   It only builds Go specific code into the Molecule object           */
00052 /*                                                                      */
00053 /************************************************************************/
00054 void Molecule::goInit() {
00055   numGoAtoms=0;
00056   energyNative=0;
00057   energyNonnative=0;
00058   atomChainTypes=NULL;
00059   goSigmaIndices=NULL;
00060   goSigmas=NULL;
00061   goWithinCutoff=NULL;
00062   goCoordinates=NULL;
00063   goResids=NULL;
00064   goPDB=NULL;
00065 }
00066 
00067 /************************************************************************/
00068 /*                                                                      */
00069 /*      FUNCTION build_go_params                                        */
00070 /*                                                                      */
00071 /*   INPUTS:                                                            */
00072 /*        fname - name of the parameter file to read                    */
00073 /*                                                                      */
00074 /*        This function reads in multiple Go parameter files            */
00075 /*   from a StringList and exhaustively processes them.                 */
00076 /*                                                                      */
00077 /************************************************************************/
00078 // JE - read in a Go parameter file
00079 void Molecule::build_go_params(StringList *g) {
00080   iout << iINFO << "Building Go Parameters" << "\n" << endi;
00081 #ifdef MEM_OPT_VERSION
00082   NAMD_die("Go forces are not supported in memory-optimized builds.");
00083 #else
00084   build_lists_by_atom();
00085 #endif
00086   int iterator = 0;
00087     do
00088     {
00089       iout << iINFO << "Reading Go File: " << iterator << "\n" << endi;
00090       read_go_file(g->data);
00091       g = g->next;
00092       iterator++;
00093     } while ( g != NULL && iterator < 100);    
00094 }
00095 
00096 // Ported by JLai from NAMD 2.7
00097 /************************************************************************/
00098 /*                                                                      */
00099 /*      FUNCTION read_go_file                                           */
00100 /*                                                                      */
00101 /*   INPUTS:                                                            */
00102 /*        fname - name of the parameter file to read                    */
00103 /*                                                                      */
00104 /*        This function reads in a Go parameter file and adds the       */ 
00105 /*   parameters from this file to the current group of parameters.      */
00106 /*   The basic logic of the routine is to first find out what type of   */
00107 /*   parameter we have in the file. Then look at each line in turn      */
00108 /*   and call the appropriate routine to add the parameters until we hit*/
00109 /*   a new type of parameter or EOF.                                    */
00110 /*                                                                      */
00111 /************************************************************************/
00112 // JE - read in a Go parameter file
00113 void Molecule::read_go_file(char *fname)
00114 
00115 {
00116 
00117   int i;                   // Counter
00118   int j;                   // Counter
00119   int  par_type=0;         //  What type of parameter are we currently
00120                            //  dealing with? (vide infra)
00121   // JLai -- uncommented
00122   int  skipline;           //  skip this line?
00123   int  skipall = 0;        //  skip rest of file;
00124   char buffer[512];           //  Buffer to store each line of the file
00125   char first_word[512];           //  First word of the current line
00126   int read_count = 0;      //  Count of input parameters on a given line
00127   int chain1 = 0;          //  First chain type for pair interaction
00128   int chain2 = 0;          //  Second chain type for pair interaction
00129   int int1;                //  First parameter int
00130   int int2;                //  Second parameter int
00131   Real r1;                 //  Parameter Real
00132   char in2[512];           //  Second parameter word
00133   GoValue *goValue1 = NULL;    //  current GoValue for loading parameters
00134   GoValue *goValue2 = NULL;    //  other current GoValue for loading parameters
00135   Bool sameGoChain = FALSE;    //  whether the current GoValue is within a chain or between chains
00136   int restrictionCount = 0;    //  number of Go restrictions set for a given chain pair
00137   FILE *pfile;                 //  File descriptor for the parameter file
00138 
00139   /*  Check to make sure that we haven't previously been told     */
00140   /*  that all the files were read                                */
00141   /*if (AllFilesRead)
00142     {
00143     NAMD_die("Tried to read another parameter file after being told that all files were read!");
00144     }*/
00145   
00146   /*  Initialize go_indices  */
00147   for (i=0; i<MAX_GO_CHAINS+1; i++) {
00148     go_indices[i] = -1;
00149   }
00150 
00151   /*  Initialize go_array   */
00152   for (i=0; i<MAX_GO_CHAINS*MAX_GO_CHAINS; i++) {
00153     go_array[i].epsilon = 1.25;
00154     go_array[i].exp_a = 12;
00155     go_array[i].exp_b = 6;
00156     go_array[i].exp_rep = 12;
00157     go_array[i].sigmaRep = 2.25;
00158     go_array[i].epsilonRep = 0.03;
00159     go_array[i].cutoff = 4.0;
00160     for (j=0; j<MAX_RESTRICTIONS; j++) {
00161       go_array[i].restrictions[j] = -1;
00162     }
00163   }
00164 
00165   /*  Try and open the file                                        */
00166   if ( (pfile = fopen(fname, "r")) == NULL)
00167     {
00168       char err_msg[256];
00169       
00170       sprintf(err_msg, "UNABLE TO OPEN GO PARAMETER FILE %s\n", fname);
00171       NAMD_die(err_msg);
00172     }
00173   
00174   /*  Keep reading in lines until we hit the EOF                        */
00175   while (NAMD_read_line(pfile, buffer) != -1)
00176     {
00177       /*  Get the first word of the line                        */
00178       NAMD_find_first_word(buffer, first_word);
00179       skipline=0;
00180       
00181       /*  First, screen out things that we ignore.                   */   
00182       /*  blank lines, lines that start with '!' or '*', lines that  */
00183       /*  start with "END".                                          */
00184       if (!NAMD_blank_string(buffer) &&
00185           (strncmp(first_word, "!", 1) != 0) &&
00186           (strncmp(first_word, "*", 1) != 0) &&
00187           (strncasecmp(first_word, "END", 3) != 0))
00188         {
00189           if ( skipall ) {
00190             iout << iWARN << "SKIPPING PART OF GO PARAMETER FILE AFTER RETURN STATEMENT\n" << endi;
00191             break;
00192           }
00193           /*  Now, determine the apropriate parameter type.   */
00194           if (strncasecmp(first_word, "chaintypes", 10)==0)
00195             {
00196               read_count=sscanf(buffer, "%s %d %d\n", first_word, &int1, &int2);
00197               if (read_count != 3) {
00198                 char err_msg[512];
00199                 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);
00200                 NAMD_die(err_msg);
00201               }
00202               chain1 = int1;
00203               chain2 = int2;
00204               if (chain1 < 1 || chain1 > MAX_GO_CHAINS ||
00205                   chain2 < 1 || chain2 > MAX_GO_CHAINS) {
00206                 char err_msg[512];
00207                 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);
00208                 NAMD_die(err_msg);
00209               }
00210               if (go_indices[chain1] == -1) {
00211                 go_indices[chain1] = NumGoChains;
00212                 NumGoChains++;
00213               }
00214               if (go_indices[chain2] == -1) {
00215                 go_indices[chain2] = NumGoChains;
00216                 NumGoChains++;
00217               }
00218               if (chain1 == chain2) {
00219                 sameGoChain = TRUE;
00220               } else {
00221                 sameGoChain = FALSE;
00222               }
00223               //goValue = &go_array[(chain1 * MAX_GO_CHAINS) + chain2];
00224               goValue1 = &go_array[(chain1*MAX_GO_CHAINS) + chain2];
00225               goValue2 = &go_array[(chain2*MAX_GO_CHAINS) + chain1];
00226 #if CODE_REDUNDANT
00227               goValue1 = &go_array[(go_indices[chain1]*MAX_GO_CHAINS) + go_indices[chain2]];
00228               goValue2 = &go_array[(go_indices[chain2]*MAX_GO_CHAINS) + go_indices[chain1]];
00229 #endif
00230               restrictionCount = 0;    //  restrictionCount applies to each chain pair separately
00231             }
00232           else if (strncasecmp(first_word, "epsilonRep", 10)==0)
00233             {
00234               read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
00235               if (read_count != 2) {}
00236               goValue1->epsilonRep = r1;
00237               if (!sameGoChain) {
00238                 goValue2->epsilonRep = r1;
00239               }
00240             }
00241           else if (strncasecmp(first_word, "epsilon", 7)==0)
00242             {
00243               // Read in epsilon
00244               read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
00245               if (read_count != 2) {}
00246               goValue1->epsilon = r1;
00247               if (!sameGoChain) {
00248                 goValue2->epsilon = r1;
00249               }
00250             }
00251           else if (strncasecmp(first_word, "exp_a", 5)==0)
00252             {
00253               read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
00254               if (read_count != 2) {}
00255               goValue1->exp_a = int1;
00256               if (!sameGoChain) {
00257                 goValue2->exp_a = int1;
00258               }
00259             }
00260           else if (strncasecmp(first_word, "exp_b", 5)==0)
00261             {
00262               read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
00263               if (read_count != 2) {}
00264               goValue1->exp_b = int1;
00265               if (!sameGoChain) {
00266                 goValue2->exp_b = int1;
00267               }
00268             }
00269           else if (strncasecmp(first_word, "exp_rep", 5)==0)
00270             {
00271               read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
00272               if (read_count != 2) {}
00273               goValue1->exp_rep = int1;
00274               if (!sameGoChain) {
00275                 goValue2->exp_rep = int1;
00276               }
00277             }
00278           else if (strncasecmp(first_word, "exp_Rep", 5)==0)
00279             {
00280               read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
00281               if (read_count != 2) {}
00282               goValue1->exp_rep = int1;
00283               if (!sameGoChain) {
00284               goValue2->exp_rep = int1;
00285               }
00286             }
00287           else if (strncasecmp(first_word, "sigmaRep", 8)==0)
00288             {
00289               read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
00290               if (read_count != 2) {}
00291               goValue1->sigmaRep = r1;
00292               if (!sameGoChain) {
00293                 goValue2->sigmaRep = r1;
00294               }
00295             }
00296           else if (strncasecmp(first_word, "cutoff", 6)==0)
00297             {
00298               read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
00299               if (read_count != 2) {}
00300               goValue1->cutoff = r1;
00301               if (!sameGoChain) {
00302                 goValue2->cutoff = r1;
00303               }
00304             }
00305           else if (strncasecmp(first_word, "restriction", 10)==0)
00306             {
00307               read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
00308               if (read_count != 2) {}
00309               if (int1 < 0) {
00310                 DebugM(3, "ERROR: residue restriction value must be nonnegative.  int1=" << int1 << "\n");
00311               }
00312               if (!sameGoChain) {
00313                 //goValue2->restrictions[restrictionCount] = int1;
00314                 DebugM(3, "ERROR: residue restrictions should not be defined between two separate chains.  chain1=" << chain1 << ", chain2=" << chain2 << "\n");
00315               }
00316               else {
00317                 goValue1->restrictions[restrictionCount] = int1;
00318               }
00319               restrictionCount++;
00320             }
00321           else if (strncasecmp(first_word, "return", 4)==0)
00322             {
00323               skipall=8;
00324               skipline=1;
00325             }        
00326           else // if (par_type == 0)
00327             {
00328               /*  This is an unknown paramter.        */
00329               /*  This is BAD                                */
00330               char err_msg[512];
00331               
00332               sprintf(err_msg, "UNKNOWN PARAMETER IN GO PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00333               NAMD_die(err_msg);
00334             }
00335         }
00336       else
00337         {
00338           skipline=1;
00339         }
00340     }
00341   
00342   /*  Close the file  */
00343   fclose(pfile);
00344   
00345   return;
00346 }
00347 /*              END OF FUNCTION read_go_file      */
00348 
00349 #if CODE_REDUNDANT
00350 /************************************************************************/
00351 /*                                                                      */
00352 /*      FUNCTION get_go_cutoff                                          */
00353 /*                                                                      */
00354 /*   INPUTS:                                                            */
00355 /*     chain1 - first chain type                                        */
00356 /*     chain2 - second chain type                                       */
00357 /*                                                                      */
00358 /*  This function gets the Go cutoff for two chain types.  The cutoff   */
00359 /*   determines whether the Go force is attractive or repulsive.        */
00360 /*                                                                      */
00361 /************************************************************************/
00362 // JE
00363 Real Molecule::get_go_cutoff(int chain1, int chain2)
00364 {
00365   return go_array[MAX_GO_CHAINS*chain1 + chain2].cutoff;
00366   #if CODE_REDUNDANT
00367   return go_array[MAX_GO_CHAINS*go_indices[chain1] + go_indices[chain2]].cutoff;
00368   #endif
00369 }
00370 /*             END OF FUNCTION get_go_cutoff      */
00371 
00372 
00373 /************************************************************************/
00374 /*                                                                      */
00375 /*      FUNCTION get_go_epsilonRep                                      */
00376 /*                                                                      */
00377 /*   INPUTS:                                                            */
00378 /*     chain1 - first chain type                                        */
00379 /*     chain2 - second chain type                                       */
00380 /*                                                                      */
00381 /*  This function gets the Go epsilonRep value for two chain types.     */
00382 /*   epsilonRep is a factor in the repulsive Go force formula.          */
00383 /*                                                                      */
00384 /************************************************************************/
00385 // JE
00386 Real Molecule::get_go_epsilonRep(int chain1, int chain2)
00387 {
00388   return go_array[MAX_GO_CHAINS*chain1 + chain2].epsilonRep;
00389   #if CODE_REDUNDANT
00390   return go_array[MAX_GO_CHAINS*go_indices[chain1] + go_indices[chain2]].epsilonRep;
00391   #endif
00392 }
00393 /*             END OF FUNCTION get_go_epsilonRep      */
00394 
00395 
00396 /************************************************************************/
00397 /*                                                                      */
00398 /*      FUNCTION get_go_sigmaRep                                        */
00399 /*                                                                      */
00400 /*   INPUTS:                                                            */
00401 /*     chain1 - first chain type                                        */
00402 /*     chain2 - second chain type                                       */
00403 /*                                                                      */
00404 /*  This function gets the Go sigmaRep value for two chain types.       */
00405 /*   sigmaRep is a factor in the repulsive Go force formula.            */
00406 /*                                                                      */
00407 /************************************************************************/
00408 // JE
00409 Real Molecule::get_go_sigmaRep(int chain1, int chain2)
00410 {
00411   return go_array[MAX_GO_CHAINS*chain1 + chain2].sigmaRep;
00412 }
00413 /*             END OF FUNCTION get_go_sigmaRep        */
00414 
00415 
00416 /************************************************************************/
00417 /*                                                                      */
00418 /*      FUNCTION get_go_epsilon                                         */
00419 /*                                                                      */
00420 /*   INPUTS:                                                            */
00421 /*     chain1 - first chain type                                        */
00422 /*     chain2 - second chain type                                       */
00423 /*                                                                      */
00424 /*  This function gets the Go epsilon value for two chain types.        */
00425 /*   epsilon is a factor in the attractive Go force formula.            */
00426 /*                                                                      */
00427 /************************************************************************/
00428 // JE
00429 Real Molecule::get_go_epsilon(int chain1, int chain2)
00430 {
00431   return go_array[MAX_GO_CHAINS*chain1 + chain2].epsilon;
00432   #if CODE_REDUNDANT
00433   return go_array[MAX_GO_CHAINS*go_indices[chain1] + go_indices[chain2]].epsilon;
00434   #endif
00435 }
00436 /*             END OF FUNCTION get_go_epsilon         */
00437 
00438 
00439 /************************************************************************/
00440 /*                                                                      */
00441 /*      FUNCTION get_go_exp_a                                           */
00442 /*                                                                      */
00443 /*   INPUTS:                                                            */
00444 /*     chain1 - first chain type                                        */
00445 /*     chain2 - second chain type                                       */
00446 /*                                                                      */
00447 /*  This function gets the Go exp_a value for two chain types.          */
00448 /*   exp_a is an exponent in the repulsive term of the attractive Go    */
00449 /*   force formula.                                                     */
00450 /*                                                                      */
00451 /************************************************************************/
00452 // JE
00453 int Molecule::get_go_exp_a(int chain1, int chain2)
00454 {
00455   return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_a;
00456   #if CODE_REDUNDANT
00457   return go_array[MAX_GO_CHAINS*go_indices[chain1] + go_indices[chain2]].exp_a;
00458   #endif
00459 }
00460 /*             END OF FUNCTION get_go_exp_a        */
00461 
00462 
00463 /************************************************************************/
00464 /*                                                                      */
00465 /*      FUNCTION get_go_exp_b                                           */
00466 /*                                                                      */
00467 /*   INPUTS:                                                            */
00468 /*     chain1 - first chain type                                        */
00469 /*     chain2 - second chain type                                       */
00470 /*                                                                      */
00471 /*  This function gets the Go exp_b value for two chain types.          */
00472 /*   exp_b is an exponent in the attractive term of the attractive Go   */
00473 /*   force formula.                                                     */
00474 /*                                                                      */
00475 /************************************************************************/
00476 // JE
00477 int Molecule::get_go_exp_b(int chain1, int chain2)
00478 {
00479   return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_b;
00480   #if CODE_REDUNDANT
00481   return go_array[MAX_GO_CHAINS*go_indices[chain1] + go_indices[chain2]].exp_b;
00482   #endif
00483 }
00484 /*             END OF FUNCTION get_go_exp_b        */
00485 
00486 
00487 /************************************************************************/
00488 /*                                                                      */
00489 /*      FUNCTION get_go_exp_rep                                         */
00490 /*                                                                      */
00491 /*   INPUTS:                                                            */
00492 /*     chain1 - first chain type                                        */
00493 /*     chain2 - second chain type                                       */
00494 /*                                                                      */
00495 /*  This function gets the Go exp_rep value for two chain types.        */
00496 /*   exp_b is an exponent in the attractive term of the attractive Go   */
00497 /*   force formula.                                                     */
00498 /*                                                                      */
00499 /************************************************************************/
00500 // JE
00501 int Molecule::get_go_exp_rep(int chain1, int chain2)
00502 {
00503   return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_rep;
00504   #if CODE_REDUNDANT
00505   return go_array[MAX_GO_CHAINS*go_indices[chain1] + go_indices[chain2]].exp_rep;
00506   #endif
00507 }
00508 /*             END OF FUNCTION get_go_exp_rep      */
00509 #endif
00510 
00511 /************************************************************************/
00512 /*                                                                      */
00513 /*      FUNCTION go_restricted                                          */
00514 /*                                                                      */
00515 /*   INPUTS:                                                            */
00516 /*     chain1 - first chain type                                        */
00517 /*     chain2 - second chain type                                       */
00518 /*     rDiff - residue ID difference to check                           */
00519 /*                                                                      */
00520 /*  This function checks whether residues with IDs rDiff apart are      */
00521 /*   restricted from Go force calculation.                              */
00522 /*                                                                      */
00523 /************************************************************************/
00524 // JE
00525 Bool Molecule::go_restricted(int chain1, int chain2, int rDiff)
00526 {
00527   int i;      //  Loop counter
00528   for(i=0; i<MAX_RESTRICTIONS; i++) {
00529     if (go_array[(MAX_GO_CHAINS*chain1) + chain2].restrictions[i]  == rDiff) {
00530       return TRUE;
00531     } else if (go_array[(MAX_GO_CHAINS*chain1) + chain2].restrictions[i] == -1) {
00532       return FALSE;
00533     }
00534   }
00535   return FALSE;
00536 }
00537 /*              END OF FUNCTION go_restricted      */
00538 
00539 // Original by JE
00540 /************************************************************************/
00541 /*                  */
00542 /*      FUNCTION print_go_params      */
00543 /*                  */
00544 /*  This is a debugging routine used to print out all the Go  */
00545 /*  parameters                */
00546 /*                  */
00547 /************************************************************************/
00548 void Molecule::print_go_params()
00549 {
00550   int i;
00551   int j;
00552   int index;
00553 
00554   DebugM(3,NumGoChains << " Go PARAMETERS 3\n" \
00555          << "*****************************************" << std::endl);
00556 
00557   for (i=0; i<NumGoChains; i++) {
00558     for (j=0; j<NumGoChains; j++) {
00559       index = (i * MAX_GO_CHAINS) + j;
00560       //  Real epsilon;    // Epsilon
00561       //  Real exp_a;      // First exponent for attractive L-J term
00562       //  Real exp_b;      // Second exponent for attractive L-J term
00563       //  Real sigmaRep;   // Sigma for repulsive term
00564       //  Real epsilonRep; // Epsilon for replusive term
00565       DebugM(3,"Go index=(" << i << "," << j << ") epsilon=" << go_array[index].epsilon \
00566              << " exp_a=" << go_array[index].exp_a << " exp_b=" << go_array[index].exp_b \
00567              << " exp_rep=" << go_array[index].exp_rep << " sigmaRep=" \
00568              << go_array[index].sigmaRep << " epsilonRep=" << go_array[index].epsilonRep \
00569              << " cutoff=" << go_array[index].cutoff << std::endl);
00570     }
00571   }
00572 
00573 }
00574 // End of port -- JLai
00575 
00576 #ifndef MEM_OPT_VERSION
00577 void Molecule::build_go_sigmas(StringList *goCoordFile, 
00578                                char *cwd)
00579 {
00580   DebugM(3,"->build_go_sigmas" << std::endl);
00581   PDB *goPDB;      //  Pointer to PDB object to use
00582   int bcol = 4;      //  Column that data is in
00583   int32 chainType = 0;      //  b value from PDB file
00584   int i;         //  Loop counter
00585   int j;         //  Loop counter
00586   int resid1;    //  Residue ID for first atom
00587   int resid2;    //  Residue ID for second atom
00588   int residDiff;     //  Difference between resid1 and resid2
00589   Real sigma;    //  Sigma calculated for a Go atom pair
00590   Real atomAtomDist;     //  Distance between two atoms
00591   Real exp_a;            //  First exponent in L-J like formula
00592   Real exp_b;            //  Second exponent in L-J like formula
00593   char filename[129];    //  Filename
00594   
00595   //JLai
00596   BigReal nativeEnergy, *native;
00597   BigReal nonnativeEnergy, *nonnative;
00598   nativeEnergy = 0;
00599   nonnativeEnergy = 0;
00600   native = &nativeEnergy;
00601   nonnative = &nonnativeEnergy;
00602 
00603   long nativeContacts = 0;
00604   long nonnativeContacts = 0;
00605 
00606   //  Get the PDB object that contains the Go coordinates.  If
00607   //  the user gave another file name, use it.  Otherwise, just use
00608   //  the PDB file that has the initial coordinates.  
00609   if (goCoordFile == NULL)
00610     {
00611       //goPDB = initial_pdb;
00612       NAMD_die("Error: goCoordFile is NULL - build_go_sigmas");
00613     }
00614   else
00615   {
00616     if (goCoordFile->next != NULL)
00617       {
00618         NAMD_die("Multiple definitions of Go atoms PDB file in configuration file");
00619       }
00620     
00621     if ( (cwd == NULL) || (goCoordFile->data[0] == '/') )
00622       {
00623         strcpy(filename, goCoordFile->data);
00624       }
00625     else
00626       {
00627         strcpy(filename, cwd);
00628         strcat(filename, goCoordFile->data);
00629       }
00630     
00631     goPDB = new PDB(filename);
00632     if ( goPDB == NULL )
00633       {
00634         NAMD_die("Memory allocation failed in Molecule::build_go_sigmas");
00635       }
00636     
00637     if (goPDB->num_atoms() != numAtoms)
00638       {
00639         NAMD_die("Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
00640       }
00641   }
00642   //  Allocate the array to hold the chain types
00643   atomChainTypes = new int32[numAtoms];
00644   //  Allocate the array to hold Go atom indices into the sigma array
00645   goSigmaIndices = new int32[numAtoms];
00646   
00647   if (atomChainTypes == NULL) {
00648     NAMD_die("memory allocation failed in Molecule::build_go_sigmas");
00649   }
00650   
00651   numGoAtoms = 0;
00652   
00653   //  Loop through all the atoms and get the Go chain types
00654   for (i=0; i<numAtoms; i++) {
00655     //  Get the chainType from the occupancy field
00656     chainType = (int32)(goPDB->atom(i))->occupancy();
00657     //  Assign the chainType value
00658     if ( chainType != 0 ) {
00659       //DebugM(3,"build_go_sigmas - atom:" << i << ", chainType:" << chainType << std::endl);
00660       atomChainTypes[i] = chainType;
00661       goSigmaIndices[i] = numGoAtoms;
00662       numGoAtoms++;
00663     }
00664     else {
00665       atomChainTypes[i] = 0;
00666       goSigmaIndices[i] = -1;
00667     }
00668     //printf("CT: %d %d %d %d\n",i,numGoAtoms,atomChainTypes[i],goSigmaIndices[i]);
00669   }
00670 
00671   // Allocate the array to hold the sigma values for each Go atom pair
00672   goSigmas = new Real[numGoAtoms*numGoAtoms];
00673   goWithinCutoff = new bool[numGoAtoms*numGoAtoms];
00674   for (i=0; i<numGoAtoms; i++) {
00675     for (j=0; j<numGoAtoms; j++) {
00676       goSigmas[i*numGoAtoms + j] = -1.0;
00677       goWithinCutoff[i*numGoAtoms + j] = false;
00678     }
00679   }
00680   //  Loop through all atom pairs and calculate sigmas
00681   DebugM(3,"    numAtoms=" << numAtoms << std::endl);
00682   for (i=0; i<numAtoms; i++) {
00683     //DebugM(3,"    i=" << i << std::endl);
00684     resid1 = (goPDB->atom(i))->residueseq();
00685     //DebugM(3,"    resid1=" << resid1 << std::endl);
00686     //if ( goSigmaIndices[i] != -1) {
00687     //  goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[i]] = 0.0;
00688     //}
00689      for (j=i+1; j<numAtoms; j++) {
00690       //DebugM(3,"    j=" << j << std::endl);
00691       resid2 = (goPDB->atom(j))->residueseq();
00692       //printf("GSIi %d %d %d\n",i,numAtoms,goSigmaIndices[i]);
00693       //printf("SAN CHECK: %d\n",goSigmaIndices[37]);
00694       //printf("GSIj %d %d %d\n",j,numAtoms,goSigmaIndices[j]);
00695       //printf("ATOMS_1to4 %d\n",!atoms_1to4(i,j));
00696       //DebugM(3,"    resid2=" << resid2 << std::endl);
00697       //  if goSigmaIndices aren't defined, don't set anything in goSigmas
00698       if ( goSigmaIndices[i] != -1 && goSigmaIndices[j] != -1 && !atoms_1to4(i,j) ) {
00699         //printf("TAKING DIFFERENCE\n");
00700         residDiff = resid2 - resid1;
00701         //printf("RESIDDIFF %d\n",residDiff);
00702         if (residDiff < 0) residDiff = -residDiff;
00703         //printf("RESIDDIFF2 %d\n",residDiff);
00704         //  if this is a Go atom pair that is not restricted
00705         //    calculate sigma
00706         //  sigmas are initially set to -1.0 if the atom pair fails go_restricted
00707         //printf("CHECKING LOOPING\n");
00708         if ( atomChainTypes[i] && atomChainTypes[j] &&
00709              !(this->go_restricted(atomChainTypes[i],atomChainTypes[j],residDiff)) ) {
00710           atomAtomDist = sqrt(pow((goPDB->atom(i))->xcoor() - (goPDB->atom(j))->xcoor(), 2.0) +
00711                               pow((goPDB->atom(i))->ycoor() - (goPDB->atom(j))->ycoor(), 2.0) +
00712                               pow((goPDB->atom(i))->zcoor() - (goPDB->atom(j))->zcoor(), 2.0));
00713           if ( atomAtomDist <= this->get_go_cutoff(atomChainTypes[i],atomChainTypes[j]) ) {
00714             exp_a = this->get_go_exp_a(atomChainTypes[i],atomChainTypes[j]);
00715             exp_b = this->get_go_exp_b(atomChainTypes[i],atomChainTypes[j]);
00716             sigma = pow(static_cast<double>(exp_b/exp_a),(1.0/(exp_a-exp_b))) * atomAtomDist;
00717             goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = sigma;
00718             goSigmas[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = sigma;
00719             goWithinCutoff[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = true;
00720             goWithinCutoff[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = true;
00721             nativeContacts++;
00722           } else {
00723             goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = 0.0;
00724             goSigmas[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = 0.0;
00725             nonnativeContacts++;
00726           }
00727         } else {
00728           goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = -1.0;
00729           goSigmas[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = -1.0;
00730         }
00731       } 
00732     }
00733   }
00734 
00735   iout << iINFO << "Number of UNIQUE    native contacts: " << nativeContacts << "\n" << endi;
00736   iout << iINFO << "Number of UNIQUE nonnative contacts: " << nonnativeContacts << "\n" << endi;
00737   
00738   //  If we had to create a PDB object, delete it now
00739   if (goCoordFile != NULL) {
00740     delete goPDB;
00741   }
00742   
00743   return;
00744 }
00745 /*      END OF FUNCTION build_go_sigmas    */
00746 
00747 void Molecule::build_go_sigmas2(StringList *goCoordFile, 
00748                                char *cwd)
00749 {
00750   DebugM(3,"->build_go_sigmas2" << std::endl);
00751   PDB *goPDB;      //  Pointer to PDB object to use
00752   int bcol = 4;      //  Column that data is in
00753   int32 chainType = 0;      //  b value from PDB file
00754   int32 residType = 0;      //  resid value from PDB file
00755   int i;         //  Loop counter
00756   int j;         //  Loop counter
00757   int resid1;    //  Residue ID for first atom
00758   int resid2;    //  Residue ID for second atom
00759   int residDiff;     //  Difference between resid1 and resid2
00760   Real sigma;    //  Sigma calculated for a Go atom pair
00761   Real atomAtomDist;     //  Distance between two atoms
00762   Real exp_a;            //  First exponent in L-J like formula
00763   Real exp_b;            //  Second exponent in L-J like formula
00764   char filename[129];    //  Filename
00765   
00766   //JLai
00767   int numLJPair = 0;
00768   long nativeContacts = 0;
00769   long nonnativeContacts = 0;
00770 
00771   //  Get the PDB object that contains the Go coordinates.  If
00772   //  the user gave another file name, use it.  Otherwise, just use
00773   //  the PDB file that has the initial coordinates.  
00774   if (goCoordFile == NULL)
00775     {
00776       //goPDB = initial_pdb;
00777       NAMD_die("Error: goCoordFile is NULL - build_go_sigmas2");
00778     }
00779   else
00780   {
00781     if (goCoordFile->next != NULL)
00782       {
00783         NAMD_die("Multiple definitions of Go atoms PDB file in configuration file");
00784       }
00785     
00786     if ( (cwd == NULL) || (goCoordFile->data[0] == '/') )
00787       {
00788         strcpy(filename, goCoordFile->data);
00789       }
00790     else
00791       {
00792         strcpy(filename, cwd);
00793         strcat(filename, goCoordFile->data);
00794       }
00795     
00796     goPDB = new PDB(filename);
00797     if ( goPDB == NULL )
00798       {
00799         NAMD_die("Memory allocation failed in Molecule::build_go_sigmas2");
00800       }
00801     
00802     if (goPDB->num_atoms() != numAtoms)
00803       {
00804         NAMD_die("Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
00805       }
00806   }
00807   //  Allocate the array to hold the chain types
00808   atomChainTypes = new int32[numAtoms];
00809   //  Allocate the array to hold Go atom indices into the sigma array
00810   goSigmaIndices = new int32[numAtoms];
00811   //  Allocate the array to hold resid 
00812   goResidIndices = new int32[numAtoms];
00813 
00814   if (atomChainTypes == NULL) {
00815     NAMD_die("memory allocation failed in Molecule::build_go_sigmas2");
00816   }
00817   
00818   numGoAtoms = 0;
00819   
00820   //  Loop through all the atoms and get the Go chain types
00821   for (i=0; i<numAtoms; i++) {
00822     //  Get the chainType from the occupancy field
00823     chainType = (int32)(goPDB->atom(i))->occupancy();
00824     //  Get the resid from the resid field
00825     residType = (int32)(goPDB->atom(i))->residueseq();
00826     //  Assign the chainType value
00827     if ( chainType != 0 ) {
00828       //DebugM(3,"build_go_sigmas2 - atom:" << i << ", chainType:" << chainType << std::endl);
00829       atomChainTypes[i] = chainType;
00830       goSigmaIndices[i] = numGoAtoms;
00831       goResidIndices[i] = residType;
00832       numGoAtoms++;
00833     }
00834     else {
00835       atomChainTypes[i] = 0;
00836       goSigmaIndices[i] = -1;
00837       goResidIndices[i] = -1;
00838     }
00839   }
00840 
00841   //Define:
00842   ResizeArray<GoPair> tmpGoPair;
00843   
00844   //  Loop through all atom pairs and calculate sigmas
00845   // Store sigmas into sorted array
00846   DebugM(3,"    numAtoms=" << numAtoms << std::endl);
00847   for (i=0; i<numAtoms; i++) {
00848     resid1 = (goPDB->atom(i))->residueseq();
00849      for (j=i+1; j<numAtoms; j++) {
00850       resid2 = (goPDB->atom(j))->residueseq();
00851       if ( goSigmaIndices[i] != -1 && goSigmaIndices[j] != -1 && !atoms_1to4(i,j) ) {
00852         residDiff = resid2 - resid1;
00853         if (residDiff < 0) residDiff = -residDiff;
00854         if ( atomChainTypes[i] && atomChainTypes[j] &&
00855              !(this->go_restricted(atomChainTypes[i],atomChainTypes[j],residDiff)) ) {
00856           atomAtomDist = sqrt(pow((goPDB->atom(i))->xcoor() - (goPDB->atom(j))->xcoor(), 2.0) +
00857                               pow((goPDB->atom(i))->ycoor() - (goPDB->atom(j))->ycoor(), 2.0) +
00858                               pow((goPDB->atom(i))->zcoor() - (goPDB->atom(j))->zcoor(), 2.0));
00859           if ( atomAtomDist <= this->get_go_cutoff(atomChainTypes[i],atomChainTypes[j]) ) {
00860             exp_a = this->get_go_exp_a(atomChainTypes[i],atomChainTypes[j]);
00861             exp_b = this->get_go_exp_b(atomChainTypes[i],atomChainTypes[j]);
00862             sigma = pow(static_cast<double>(exp_b/exp_a),(1.0/(exp_a-exp_b))) * atomAtomDist;
00863             double tmpA = pow(sigma,exp_a);
00864             double tmpB = pow(sigma,exp_b);
00865             GoPair gp;
00866             GoPair gp2;
00867             gp.goIndxA = i;
00868             gp.goIndxB = j;
00869             gp.A = tmpA;
00870             gp.B = tmpB;
00871             tmpGoPair.add(gp);
00872             gp2.goIndxA = j;
00873             gp2.goIndxB = i;
00874             gp2.A = tmpA;
00875             gp2.B = tmpB;
00876             tmpGoPair.add(gp2);
00877             nativeContacts++;
00878           } else {
00879             nonnativeContacts++;
00880           }
00881         } 
00882       } 
00883     }
00884   }
00885 
00886   iout << iINFO << "Number of UNIQUE    native contacts: " << nativeContacts << "\n" << endi;
00887   iout << iINFO << "Number of UNIQUE nonnative contacts: " << nonnativeContacts << "\n" << endi;
00888 
00889   // Copies the resizeArray into a block of continuous memory
00890   std::sort(tmpGoPair.begin(),tmpGoPair.end(),goPairCompare);
00891   goNumLJPair = 2*nativeContacts;
00892   goIndxLJA = new int[goNumLJPair];
00893   goIndxLJB = new int[goNumLJPair];
00894   goSigmaPairA = new double[goNumLJPair];
00895   goSigmaPairB = new double[goNumLJPair];
00896   for(i=0; i< goNumLJPair; i++) {
00897     goIndxLJA[i] = tmpGoPair[i].goIndxA;
00898     goIndxLJB[i] = tmpGoPair[i].goIndxB;
00899     goSigmaPairA[i] = tmpGoPair[i].A;
00900     goSigmaPairB[i] = tmpGoPair[i].B;
00901   }
00902 
00903   pointerToGoBeg = new int[numAtoms];
00904   pointerToGoEnd = new int[numAtoms];
00905   int oldIndex = -1;
00906   for(i=0; i<numAtoms; i++) {
00907     pointerToGoBeg[i] = -1;
00908     pointerToGoEnd[i] = -2;
00909   }
00910   for(i=0; i < goNumLJPair; i++) {
00911     if(pointerToGoBeg[goIndxLJA[i]] == -1) {
00912       pointerToGoBeg[goIndxLJA[i]] = i;
00913       oldIndex = goIndxLJA[i];
00914     }
00915     pointerToGoEnd[oldIndex] = i;
00916   }
00917 
00918   //  If we had to create a PDB object, delete it now
00919   if (goCoordFile != NULL) {
00920     delete goPDB;
00921   }
00922     
00923   return;
00924 }
00925 /*      END OF FUNCTION build_go_sigmas2    */
00926 
00927 bool Molecule::goPairCompare(GoPair first, GoPair second) {
00928   if(first.goIndxA < second.goIndxA) {
00929     return true;
00930   } else if(first.goIndxA == second.goIndxA) {
00931     return (first.goIndxB == second.goIndxB);
00932   } 
00933   return false;
00934 }
00935 
00936     /************************************************************************/
00937     /*                                                                      */
00938     /*      JE - FUNCTION build_go_arrays                                   */
00939     /*                                                                      */
00940     /*   INPUTS:                                                            */
00941     /*  goCoordFile - Value of Go coordinate file from config file          */
00942     /*  cwd - Current working directory                                     */
00943     /*                                                                      */
00944     /*  This function builds arrays that support sigma calculation for L-J  */
00945     /* style Go forces.  It takes the name of the PDB file.  It then builds */
00946     /* an array identifying atoms to which Go forces apply.                 */
00947     /*                                                                      */
00948     /************************************************************************/
00949 // JE
00950 void Molecule::build_go_arrays(StringList *goCoordFile, 
00951                               char *cwd)
00952 {
00953   DebugM(3,"->build_go_arrays" << std::endl);
00954   //PDB *goPDB;      //  Pointer to PDB object to use
00955   int bcol = 4;      //  Column that data is in
00956   int32 chainType = 0;      //  b value from PDB file
00957   int i;         //  Loop counter
00958   int j;         //  Loop counter
00959   BigReal atomAtomDist;     //  Distance between two atoms -- JLai put back 
00960   int resid1;    //  Residue ID for first atom
00961   int resid2;    //  Residue ID for second atom
00962   int residDiff;     //  Difference between resid1 and resid2
00963   int goIndex;       //  Index into the goCoordinates array
00964   int goIndx;        //  Index into the goCoordinates array
00965   char filename[129];    //  Filename
00966   
00967   //JLai
00968   BigReal nativeEnergy, *native;
00969   BigReal nonnativeEnergy, *nonnative;
00970   nativeEnergy = 0;
00971   nonnativeEnergy = 0;
00972   native = &nativeEnergy;
00973   nonnative = &nonnativeEnergy;
00974 
00975   long nativeContacts = 0;
00976   long nonnativeContacts = 0;
00977 
00978   //  Get the PDB object that contains the Go coordinates.  If
00979   //  the user gave another file name, use it.  Otherwise, just use
00980   //  the PDB file that has the initial coordinates.  
00981   if (goCoordFile == NULL)
00982     {
00983       //goPDB = initial_pdb;
00984       NAMD_die("Error: goCoordFile is NULL - build_go_arrays");
00985     }
00986   else
00987   {
00988     if (goCoordFile->next != NULL)
00989       {
00990         NAMD_die("Multiple definitions of Go atoms PDB file in configuration file");
00991       }
00992     
00993     if ( (cwd == NULL) || (goCoordFile->data[0] == '/') )
00994       {
00995         strcpy(filename, goCoordFile->data);
00996       }
00997     else
00998       {
00999         strcpy(filename, cwd);
01000         strcat(filename, goCoordFile->data);
01001       }
01002     
01003     goPDB = new PDB(filename);
01004     if ( goPDB == NULL )
01005       {
01006         NAMD_die("goPDB memory allocation failed in Molecule::build_go_arrays");
01007       }
01008     
01009     if (goPDB->num_atoms() != numAtoms)
01010       {
01011         NAMD_die("Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
01012       }
01013   }
01014   
01015   //  Allocate the array to hold Go atom indices into the sigma array
01016   goSigmaIndices = new int32[numAtoms];
01017   
01018   if (goSigmaIndices == NULL) {
01019     NAMD_die("goSigmaIndices memory allocation failed in Molecule::build_go_arrays");
01020   }
01021   
01022   numGoAtoms = 0;
01023   
01024   //  Loop through all the atoms and get the Go chain types
01025   for (i=0; i<numAtoms; i++) {
01026     chainType = (int32)(goPDB->atom(i))->occupancy();
01027     if ( chainType != 0 ) {
01028       DebugM(3,"build_go_arrays - atom:" << i << std::endl);
01029       goSigmaIndices[i] = numGoAtoms;
01030       numGoAtoms++;
01031     }
01032     else {
01033       goSigmaIndices[i] = -1;
01034     }
01035   }
01036 
01037   // Allocate the array to hold the sigma values for each Go atom pair
01049   //  Allocate the array to hold the chain types
01050   atomChainTypes = new int32[numGoAtoms];
01051 
01052   if (atomChainTypes == NULL) {
01053     NAMD_die("atomChainTypes memory allocation failed in Molecule::build_go_arrays");
01054   }
01055 
01056   // Allocate the array to hold (x,y,z) coordinates for all of the Go atoms
01057   goCoordinates = new Real[numGoAtoms*3];
01058 
01059   if (goCoordinates == NULL) {
01060     NAMD_die("goCoordinates memory allocation failed in Molecule::build_go_arrays");
01061   }
01062 
01063   goResids = new int[numGoAtoms];
01064 
01065   // Allocate the array to hold PDB residu IDs for all of the Go atoms
01066   if (goResids == NULL) {
01067     NAMD_die("goResids memory allocation failed in Molecule::build_go_arrays");
01068   }
01069   
01070   for (i=0; i<numAtoms; i++) {
01071     goIndex = goSigmaIndices[i];
01072     if (goIndex != -1) {
01073       //  Assign the chainType value!
01074       //  Get the chainType from the occupancy field
01075       atomChainTypes[goIndex] = (int32)(goPDB->atom(i))->occupancy();
01076       goCoordinates[goIndex*3] = goPDB->atom(i)->xcoor();
01077       goCoordinates[goIndex*3 + 1] = goPDB->atom(i)->ycoor();
01078       goCoordinates[goIndex*3 + 2] = goPDB->atom(i)->zcoor();
01079       goResids[goIndex] = goPDB->atom(i)->residueseq();
01080     }
01081   }
01082       // JLai
01083   energyNative = 0;
01084   energyNonnative = 0;
01085   //printf("INIT ENERGY: (N) %f (NN) %f\n", energyNative, energyNonnative);
01086   for (i=0; i<numAtoms-1; i++) {
01087     goIndex = goSigmaIndices[i];
01088     if (goIndex != -1) {
01089       for (j=i+1; j<numAtoms;j++) {
01090         goIndx = goSigmaIndices[j];
01091         if (goIndx != -1) {
01092           resid1 = (goPDB->atom(i))->residueseq();
01093           resid2 = (goPDB->atom(j))->residueseq();
01094           residDiff = resid2 - resid1;
01095           if (residDiff < 0) residDiff = -residDiff;
01096           if (atomChainTypes[goIndex] && atomChainTypes[goIndx] &&
01097               !(this->go_restricted(atomChainTypes[goIndex],atomChainTypes[goIndx],residDiff)) &&
01098               !atoms_1to4(i,j)) {
01099             atomAtomDist = sqrt(pow((goPDB->atom(i))->xcoor() - (goPDB->atom(j))->xcoor(), 2.0) +
01100                                 pow((goPDB->atom(i))->ycoor() - (goPDB->atom(j))->ycoor(), 2.0) +
01101                                 pow((goPDB->atom(i))->zcoor() - (goPDB->atom(j))->zcoor(), 2.0));
01102             if (atomAtomDist <= this->get_go_cutoff(atomChainTypes[goIndex],atomChainTypes[goIndx]) ) {
01103               nativeContacts++;
01104             } else {
01105               nonnativeContacts++;
01106             }    
01107           }
01108         }
01109       }
01110     }
01111   }
01112   iout << iINFO << "Number of UNIQUE    native contacts: " << nativeContacts     << "\n" << endi;
01113   iout << iINFO << "Number of UNIQUE nonnative contacts: " << nonnativeContacts  << "\n" << endi;
01114   
01115   //  If we had to create a PDB object, delete it now
01116   if (goCoordFile != NULL) {
01117     delete goPDB;
01118   }
01119   
01120   return;
01121 }
01122 /*      END OF FUNCTION build_go_arrays    */
01123 #endif // #ifndef MEM_OPT_VERSION
01124 
01125 /************************************************************************/
01126 /*                                                                      */
01127 /*      FUNCTION print_go_sigmas                                        */
01128 /*                                                                      */
01129 /*  print_go_sigmas prints out goSigmas, the array of sigma parameters  */
01130 /*   used in the L-J type Go force calculations                         */
01131 /*                                                                      */
01132 /************************************************************************/
01133 // JE
01134 void Molecule::print_go_sigmas()
01135 {
01136   int i;  //  Counter
01137   int j;  //  Counter
01138   Real sigma;
01139 
01140   DebugM(3,"GO SIGMA ARRAY\n" \
01141          << "***************************" << std::endl);
01142   DebugM(3, "numGoAtoms: " << numGoAtoms << std::endl);
01143 
01144   if (goSigmaIndices == NULL) {
01145     DebugM(3, "GO SIGMAS HAVE NOT BEEN BUILT" << std::endl);
01146     return;
01147   }
01148 
01149   for (i=0; i<numAtoms; i++) {
01150     for (j=0; j<numAtoms; j++) {
01151       if ( goSigmaIndices[i] != -1 && goSigmaIndices[j] != -1 ) {
01152         //DebugM(3, "i: " << i << ", j: " << j << std::endl);
01153         sigma = goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]];
01154         if (sigma > 0.0) {
01155           DebugM(3, "(" << i << "," << j << ") - +" << sigma << " ");
01156         }
01157         else {
01158           DebugM(3, "(" << i << "," << j << ") - " << sigma << " ");
01159         }
01160       } else {
01161         //DebugM(3, "0 ");
01162       }
01163     }
01164     if ( goSigmaIndices[i] != -1 ) {
01165       DebugM(3, "-----------" << std::endl);
01166     }
01167   }
01168   return;
01169 }
01170 /*      END OF FUNCTION print_go_sigmas     */
01171 
01172 // JLai
01173 BigReal Molecule::get_gro_force2(BigReal x,
01174                                  BigReal y,
01175                                  BigReal z,
01176                                  int atom1,
01177                                  int atom2,
01178                                  BigReal* pairLJEnergy,
01179                                  BigReal* pairGaussEnergy) const
01180 {
01181   //Initialize return energies to zero
01182   *pairLJEnergy = 0.0;
01183   *pairGaussEnergy = 0.0;
01184 
01185   // Linear search for LJ data
01186   int LJIndex = -1;
01187   int LJbegin = pointerToLJBeg[atom1];
01188   int LJend = pointerToLJEnd[atom1];
01189   for(int i = LJbegin; i <= LJend; i++) {
01190     if(indxLJB[i] == atom2) {
01191       LJIndex = i;
01192       break;
01193     }
01194   }
01195 
01196   // Linear search for Gaussian data
01197   int GaussIndex = -1;
01198   int Gaussbegin = pointerToGaussBeg[atom1];
01199   int Gaussend = pointerToGaussEnd[atom1];
01200   for(int i = Gaussbegin; i <= Gaussend; i++) {
01201     if(indxGaussB[i] == atom2) {
01202       GaussIndex = i;
01203       break;
01204     }
01205   }
01206 
01207   if( LJIndex == -1 && GaussIndex == -1) {
01208     return 0;
01209   } else {
01210     // Code to calculate distances because the pair was found in one of the lists
01211     BigReal r2 = x*x + y*y + z*z;
01212     BigReal r = sqrt(r2);
01213     BigReal ri = 1/r;
01214     BigReal ri6 = (ri*ri*ri*ri*ri*ri);
01215     BigReal ri12 = ri6*ri6;
01216     BigReal ri13 = ri12*ri;
01217     BigReal LJ = 0;
01218     BigReal Gauss = 0;
01219     // Code to calculate LJ 
01220     if (LJIndex != -1) {
01221       BigReal ri7 = ri6*ri;
01222       LJ = (12*(pairC12[LJIndex]*ri13) - 6*(pairC6[LJIndex]*ri7));
01223       *pairLJEnergy = (pairC12[LJIndex]*ri12 - pairC6[LJIndex]*ri6);
01224       //std::cout << pairC12[LJIndex] << " " << pairC6[LJIndex] << " " << ri13 << " " << ri7 << " " << LJ << " " << r << "\n";
01225     }
01226     // Code to calculate Gaussian
01227     if (GaussIndex != -1) {
01228       BigReal gr = 12*gRepulsive[GaussIndex]*ri13;
01229       BigReal r1prime = r - gMu1[GaussIndex];
01230       BigReal tmp1 = r1prime * r1prime;
01231       BigReal r2prime = r - gMu2[GaussIndex];
01232       BigReal tmp2 = r2prime * r2prime;
01233       BigReal tmp_gauss1 = 0;
01234       BigReal one_gauss1 = 1;
01235       BigReal tmp_gauss2 = 0;
01236       BigReal one_gauss2 = 1;
01237       if (giSigma1[GaussIndex] != 0) {
01238         tmp_gauss1 = exp(-tmp1*giSigma1[GaussIndex]);
01239         one_gauss1 = 1 - tmp_gauss1;
01240       }
01241       if (giSigma2[GaussIndex] != 0) {
01242         tmp_gauss2 = exp(-tmp2*giSigma2[GaussIndex]);
01243         one_gauss2 = 1 - tmp_gauss2;
01244       } 
01245       BigReal A = gA[GaussIndex];
01246       Gauss = gr*one_gauss1*one_gauss2 - 2*A*tmp_gauss1*one_gauss2*r1prime*giSigma1[GaussIndex] \
01247         - 2*tmp_gauss1*one_gauss2*r1prime*giSigma1[GaussIndex]*gRepulsive[GaussIndex]*ri12 - 2*A*tmp_gauss2*one_gauss1*r2prime*giSigma2[GaussIndex] \
01248         - 2*tmp_gauss2*one_gauss1*r2prime*giSigma2[GaussIndex]*gRepulsive[GaussIndex]*ri12;
01249       *pairGaussEnergy = A*(-1+(one_gauss1)*(one_gauss2)*(1+gRepulsive[GaussIndex]*ri12/A));
01250     }
01251     //std::cout << "Net force: " << (LJ + Gauss) << " with ri " << (LJ + Gauss)*ri << "\n";
01252     return (LJ + Gauss)*ri;
01253   }
01254   return 0;
01255 }
01256 // End of get_gro_force2
01257 // JLai
01258 
01259 // JE
01260 BigReal Molecule::get_go_force(BigReal r, 
01261                                int atom1,
01262                                int atom2,
01263                                BigReal* goNative,
01264                                BigReal* goNonnative) const
01265 {
01266   BigReal goForce = 0.0;
01267   Real pow1;
01268   Real pow2;
01269   //  determine which Go chain pair we are working with
01270   //DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ")" << std::endl);
01271   int32 chain1 = atomChainTypes[atom1];
01272   int32 chain2 = atomChainTypes[atom2];
01273 
01274   //DebugM(3,"  chain1:" << chain1 << ", chain2:" << chain2 << std::endl);
01275   if (chain1 == 0 || chain2 == 0)  return 0.0;
01276 
01277   //  retrieve Go cutoff for this chain pair
01278   //TMP// JLai -- I'm going to replace this with a constant accessor.  This is just a temporary thing
01279   Real goCutoff = const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2);
01280   //DebugM(3,"  goCutoff:" << goCutoff << std::endl);
01281   if (goCutoff == 0)  return 0.0;
01282   //  if repulsive then calculate repulsive
01283   //  sigmas are initially set to -1.0 if the atom pair fails go_restricted
01284   if (goSigmas[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]] != -1.0) {
01285     if (!goWithinCutoff[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]]) {
01286       Real epsilonRep = const_cast<Molecule*>(this)->get_go_epsilonRep(chain1,chain2);
01287       Real sigmaRep = const_cast<Molecule*>(this)->get_go_sigmaRep(chain1,chain2);
01288       int exp_rep = const_cast<Molecule*>(this)->get_go_exp_rep(chain1,chain2);
01289       pow1 = pow(sigmaRep/r,exp_rep);
01290       goForce = 4*((exp_rep/(r*r)) * epsilonRep * pow1);
01291       *goNative = 0.0;
01292       *goNonnative = (4 * epsilonRep * pow1 );
01293       //DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", epsilonRep:" << epsilonRep << ", sigmaRep:" << sigmaRep << ", r:" << r << ", goForce:" << goForce << std::endl);
01294     }
01295     //  if attractive then calculate attractive
01296     else {
01297       int goSigmaIndex1 = goSigmaIndices[atom1];
01298       int goSigmaIndex2 = goSigmaIndices[atom2];
01299       if (goSigmaIndex1 != -1 && goSigmaIndex2 != -1) {
01300         Real epsilon = const_cast<Molecule*>(this)->get_go_epsilon(chain1,chain2);
01301         int exp_a = const_cast<Molecule*>(this)->get_go_exp_a(chain1,chain2);
01302         int exp_b = const_cast<Molecule*>(this)->get_go_exp_b(chain1,chain2);
01303         Real sigma_ij = goSigmas[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]];
01304         // Positive gradient of potential, not negative gradient of potential
01305         pow1 = pow(sigma_ij/r,exp_a);
01306         pow2 = pow(sigma_ij/r,exp_b);
01307         goForce = ((4/(r*r)) * epsilon * (exp_a * pow1 - exp_b * pow2));
01308         //DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", sigma_ij:" << sigma_ij << ", r:" << r << ", goForce:" << goForce << std::endl);
01309         *goNative = (4 * epsilon * ( pow1 -  pow2 ) );
01310         *goNonnative = 0.0;
01311       }
01312     }
01313   }
01314   //DebugM(3,"goForce:" << goForce << std::endl);
01315   return goForce;
01316 }
01317 /*      END OF FUNCTION get_go_force_old       */
01318 
01319 
01320     /************************************************************************/
01321     /*                                                                      */
01322     /*      JE - FUNCTION get_go_force_new                                  */
01323     /*                                                                      */
01324     /*   INPUTS:                                                            */
01325     /*  r - distance between the two atoms                                  */
01326     /*  atom1 - the ID of the first atom                                    */
01327     /*  atom2 - the ID of the second atom                                   */
01328     /*                                                                      */
01329     /*  This function calculates the Go force between two atoms.  If the    */
01330     /*   atoms do not have Go parameters or sigmas, 0 is returned.          */
01331     /*                                                                      */
01332     /************************************************************************/
01333 // JE
01334 BigReal Molecule::get_go_force_new(BigReal r,
01335                                    int atom1,
01336                                    int atom2,
01337                                    BigReal* goNative,
01338                                    BigReal* goNonnative) const
01339 {
01340   int resid1;
01341   int resid2;
01342   int residDiff;
01343   Real xcoorDiff;
01344   Real ycoorDiff;
01345   Real zcoorDiff;
01346   Real atomAtomDist;
01347   Real exp_a;
01348   Real exp_b;
01349   Real sigma_ij;
01350   Real epsilon;
01351   Real epsilonRep;
01352   Real sigmaRep;
01353   Real expRep;
01354   Real pow1;
01355   Real pow2;
01356   
01357   BigReal goForce = 0.0;
01358   *goNative = 0;
01359   *goNonnative = 0;
01360 
01361   //  determine which Go chain pair we are working with
01362   DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ")" << std::endl);
01363   int goIndex1 = goSigmaIndices[atom1];
01364   int goIndex2 = goSigmaIndices[atom2];
01365 
01366   int32 chain1 = atomChainTypes[goIndex1];
01367   int32 chain2 = atomChainTypes[goIndex2];
01368 
01369   DebugM(3,"  chain1:" << chain1 << ", chain2:" << chain2 << std::endl);
01370   if (chain1 == 0 || chain2 == 0)  return 0.0;
01371 
01372   //  retrieve Go cutoff for this chain pair
01373   Real goCutoff = const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2);
01374   DebugM(3,"  goCutoff:" << goCutoff << std::endl);
01375   if (goCutoff == 0)  return 0.0;
01376 
01377   //  sigmas are initially set to -1.0 if the atom pair fails go_restricted
01378   //  no goSigmas array anymore
01379   //Real sigma_ij = goSigmas[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]];
01380 
01381   // XXX - used to be a condition for the following if
01382   //if the atoms are within 4 of each other
01383   //if ( !atoms_1to4(atom1,atom2) ) {
01384 
01385   //  if goSigmaIndices aren't defined, don't calculate forces
01386   if ( goIndex1 != -1 && goIndex2 != -1 ) {
01387     resid1 = goResids[goIndex1];
01388     resid2 = goResids[goIndex2];
01389     residDiff = resid2 - resid1;
01390     if (residDiff < 0) residDiff = -residDiff;
01391     //  if this is a Go atom pair that is not restricted
01392     if ( !(const_cast<Molecule*>(this)->go_restricted(chain1,chain2,residDiff)) ) {
01393       xcoorDiff = goCoordinates[goIndex1*3] - goCoordinates[goIndex2*3];
01394       ycoorDiff = goCoordinates[goIndex1*3 + 1] - goCoordinates[goIndex2*3 + 1];
01395       zcoorDiff = goCoordinates[goIndex1*3 + 2] - goCoordinates[goIndex2*3 + 2];
01396       atomAtomDist = sqrt(xcoorDiff*xcoorDiff + ycoorDiff*ycoorDiff + zcoorDiff*zcoorDiff);
01397       
01398       //  if attractive then calculate attractive
01399       if ( atomAtomDist <= const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2) ) {
01400         exp_a = const_cast<Molecule*>(this)->get_go_exp_a(chain1,chain2);
01401         exp_b = const_cast<Molecule*>(this)->get_go_exp_b(chain1,chain2);
01402         sigma_ij = pow(static_cast<double>(exp_b/exp_a),(1.0/(exp_a-exp_b))) * atomAtomDist;
01403         
01404         // [JLai] print out atoms involved in native contacts
01405         // printf("ATOM1: %d C1: %d ATOM2: %d C2: %d\n", atom1,chain1,atom2,chain2);
01406 
01407         epsilon = const_cast<Molecule*>(this)->get_go_epsilon(chain1,chain2);
01408         pow1 = pow(sigma_ij/r,static_cast<double>(exp_a));
01409         pow2 = pow(sigma_ij/r,static_cast<double>(exp_b));
01410         //goForce = ((4/r) * epsilon * (exp_a * pow(sigma_ij/r,exp_a) - exp_b * pow(sigma_ij/r,exp_b)));
01411         goForce = ((4/(r*r)) * epsilon * (exp_a * pow1 - exp_b * pow2));
01412         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);
01413         //goEnergy = (4 * epsilon * ( pow(sigma_ij/r,exp_a) -  pow(sigma_ij/r,exp_b) ) ); // JLai I changed some of the expressions
01414         *goNative = (4 * epsilon * ( pow1 -  pow2 ) ); 
01415         *goNonnative = 0;
01416       }
01417       
01418       //  if repulsive then calculate repulsive
01419       else {
01420         epsilonRep = const_cast<Molecule*>(this)->get_go_epsilonRep(chain1,chain2);
01421         sigmaRep = const_cast<Molecule*>(this)->get_go_sigmaRep(chain1,chain2);
01422         expRep = const_cast<Molecule*>(this)->get_go_exp_rep(chain1,chain2);
01423         pow1 = pow(sigmaRep/r,(BigReal)expRep);
01424         //goForce = ((12.0/r) * epsilonRep * pow(sigmaRep/r,12.0));
01425         goForce = (4*(expRep/(r*r)) * epsilonRep * pow1);
01426         DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", epsilonRep:" << epsilonRep << ", sigmaRep:" << sigmaRep << ", r:" << r << ", goForce:" << goForce << std::endl);
01427         //goEnergy = (4 * epsilonRep * pow(sigmaRep/r,12.0)); // JLai I changed some of the expressions
01428         *goNonnative = (4 * epsilonRep * pow1); 
01429         *goNative = 0;
01430       }
01431     }
01432   }
01433   
01434   //DebugM(3,"goForce:" << goForce << std::endl);
01435   return goForce;
01436 }
01437 /*      END OF FUNCTION get_go_force_new   */
01438 
01439 
01440     /************************************************************************/
01441     /*                                                                      */
01442     /*   JLai - FUNCTION get_go_force2                                      */
01443     /*                                                                      */
01444     /*   INPUTS:                                                            */
01445     /*  x - the x distance between p_i and p_j                              */
01446     /*  y - the y distance between p_i and p_j                              */
01447     /*  z - the z distance between p_i and p_j                              */
01448     /*  atom1 - the ID of the second atom                                   */
01449     /*  atom2 - the ID of the second atom                                   */
01450     /*                                                                      */
01451     /*  This function returns the force between two input atoms give their  */
01452 /*  distance and their atom indices.                                        */
01453     /*                                                                      */
01454     /************************************************************************/
01455 // JLai
01456 BigReal Molecule::get_go_force2(BigReal x,
01457                                 BigReal y,
01458                                 BigReal z,
01459                                 int atom1,
01460                                 int atom2,
01461                                 BigReal *goNative,
01462                                 BigReal *goNonnative) const
01463 {
01464  
01465   // Check to see if restricted.  If so, escape function early
01466   int32 chain1 = atomChainTypes[atom1];
01467   int32 chain2 = atomChainTypes[atom2];
01468 
01469   if(chain1 == 0 || chain2 == 0) return 0.0;
01470   Molecule *mol = const_cast<Molecule*>(this);
01471   Real goCutoff = mol->get_go_cutoff(chain1,chain2);
01472   if(goCutoff == 0) return 0.0;
01473 
01474   int resid1 = goResidIndices[atom1];
01475   int resid2 = goResidIndices[atom2];
01476   int residDiff = abs(resid1 - resid2);
01477   if((mol->go_restricted(chain1,chain2,residDiff))) {
01478     return 0.0;
01479   }
01480 
01481   int LJIndex = -1;
01482   int LJbegin = pointerToGoBeg[atom1];
01483   int LJend = pointerToGoEnd[atom1];
01484   for(int i = LJbegin; i <= LJend; i++) {
01485     if(goIndxLJB[i] == atom2) {
01486       LJIndex = i;
01487     }
01488   }
01489   
01490   BigReal r2 = x*x + y*y + z*z;
01491   BigReal r = sqrt(r2);
01492 
01493   if (LJIndex == -1) {
01494     int exp_rep = const_cast<Molecule*>(this)->get_go_exp_rep(chain1,chain2);
01495     BigReal epsilonRep = const_cast<Molecule*>(this)->get_go_epsilonRep(chain1, chain2);
01496     BigReal sigmaRep = const_cast<Molecule*>(this)->get_go_sigmaRep(chain1, chain2);
01497     double sigmaRepPow = pow(sigmaRep,exp_rep);
01498     BigReal LJ = (4*epsilonRep*exp_rep*sigmaRepPow*pow(r,-(exp_rep+1)));
01499     *goNative = 0;
01500     *goNonnative = (4*epsilonRep*sigmaRepPow*pow(r,-exp_rep));
01501     //*goNonnative = (4*epsilonRep * pow(sigmaRep/r,exp_rep));
01502     return (LJ/r);
01503   } else {
01504     // Code to calculate distances because the pair was found in one of the lists
01505     int exp_a = const_cast<Molecule*>(this)->get_go_exp_a(chain1,chain2);
01506     int exp_b = const_cast<Molecule*>(this)->get_go_exp_b(chain1,chain2);
01507     // We want the force, so we have to take the n+1 derivative
01508     BigReal powA = pow(r,-(exp_a + 1));
01509     BigReal powB = pow(r,-(exp_b + 1));
01510     BigReal powaa = pow(r,-exp_a);
01511     BigReal powbb = pow(r,-exp_b);
01512     BigReal epsilon = const_cast<Molecule*>(this)->get_go_epsilon(chain1,chain2);
01513     BigReal LJ = 4 * epsilon * (exp_a*goSigmaPairA[LJIndex]*powA - exp_b*goSigmaPairB[LJIndex]*powB);
01514     *goNative =  4 * epsilon * (goSigmaPairA[LJIndex]*powaa - goSigmaPairB[LJIndex]*powbb);
01515     *goNonnative = 0;
01516     return (LJ/r);
01517   }
01518   return 0;
01519 }
01520 // JLai
01521 /*      END OF FUNCTION get_go_force2   */
01522 
01523 #ifndef MEM_OPT_VERSION
01524     /************************************************************************/
01525     /*                                                                      */
01526     /*      JE - FUNCTION atoms_1to4                                        */
01527     /*                                                                      */
01528     /*   INPUTS:                                                            */
01529     /*  atom1 - the ID of the first atom                                  */
01530     /*  atom2 - the ID of the second atom                                 */
01531     /*                                                                      */
01532     /*  This function tells whether or not the two input atoms are within   */
01533     /*   1 to 4 bonds away from each other: bonds, angles, or dihedrals.    */
01534     /*                                                                      */
01535     /************************************************************************/
01536 // JE
01537 Bool Molecule::atoms_1to4(unsigned int atom1,
01538                           unsigned int atom2)
01539 {
01540   int bondNum;   //  Bonds in bonded list
01541   int angleNum;  //  Angles in angle list
01542   int dihedralNum;   //  Dihedrals in dihedral list
01543   int *bonds;
01544   int *angles;
01545   int *dihedrals;
01546   Bond *bond;     //  Temporary bond structure
01547   Angle *angle;   //  Temporary angle structure
01548   Dihedral *dihedral; //  Temporary dihedral structure
01549 
01550   DebugM(2,"atoms_1to4(" << atom1 << "," << atom2 << ")" << std::endl);
01551 
01552   bonds = get_bonds_for_atom(atom1);
01553   bondNum = *bonds;
01554   while(bondNum != -1) {
01555     bond = get_bond(bondNum);
01556     DebugM(2,"bond  atom1:" << bond->atom1 << ", atom2:" << bond->atom2 << std::endl);
01557     if (atom2 == bond->atom1 || atom2 == bond->atom2) {
01558       return TRUE;
01559     }
01560     bondNum = *(++bonds);
01561   }
01562 
01563   bonds = get_bonds_for_atom(atom2);
01564   bondNum = *bonds;
01565   while(bondNum != -1) {
01566     bond = get_bond(bondNum);
01567     DebugM(2,"bond  atom1:" << bond->atom1 << ", atom2:" << bond->atom2 << std::endl);
01568     if (atom1 == bond->atom1 || atom1 == bond->atom2) {
01569       return TRUE;
01570     }
01571     bondNum = *(++bonds);
01572   }
01573 
01574   angles = get_angles_for_atom(atom1);
01575   angleNum = *angles;
01576   while(angleNum != -1) {
01577     angle = get_angle(angleNum);
01578     DebugM(2,"angle  atom1:" << angle->atom1 << ", atom2:" << angle->atom2 << ", atom3:" << angle->atom3 << std::endl);
01579     if (atom2 == angle->atom1 || atom2 == angle->atom2 || atom2 == angle->atom3) {
01580       return TRUE;
01581     }
01582     angleNum = *(++angles);
01583   }
01584 
01585   angles = get_angles_for_atom(atom2);
01586   angleNum = *angles;
01587   while(angleNum != -1) {
01588     angle = get_angle(angleNum);
01589     DebugM(2,"angle  atom1:" << angle->atom1 << ", atom2:" << angle->atom2 << ", atom3:" << angle->atom3 << std::endl);
01590     if (atom1 == angle->atom1 || atom1 == angle->atom2 || atom1 == angle->atom3) {
01591       return TRUE;
01592     }
01593     angleNum = *(++angles);
01594   }
01595 
01596   dihedrals = get_dihedrals_for_atom(atom1);
01597   dihedralNum = *dihedrals;
01598   while(dihedralNum != -1) {
01599     dihedral = get_dihedral(dihedralNum);
01600     DebugM(2,"dihedral  atom1:" << dihedral->atom1 << ", atom2:" << dihedral->atom2 << ", atom3:" << dihedral->atom3 << ", atom4:" << dihedral->atom4 << std::endl);
01601     if (atom2 == dihedral->atom1 || atom2 == dihedral->atom2 \
01602         || atom2 == dihedral->atom3 || atom2 == dihedral->atom4) {
01603       return TRUE;
01604     }
01605     dihedralNum = *(++dihedrals);
01606   }
01607 
01608   dihedrals = get_dihedrals_for_atom(atom2);
01609   dihedralNum = *dihedrals;
01610   while(dihedralNum != -1) {
01611     dihedral = get_dihedral(dihedralNum);
01612     DebugM(2,"dihedral  atom1:" << dihedral->atom1 << ", atom2:" << dihedral->atom2 << ", atom3:" << dihedral->atom3 << ", atom4:" << dihedral->atom4 << std::endl);
01613     if (atom1 == dihedral->atom1 || atom1 == dihedral->atom2 \
01614         || atom1 == dihedral->atom3 || atom1 == dihedral->atom4) {
01615       return TRUE;
01616     }
01617     dihedralNum = *(++dihedrals);
01618   }
01619   
01620   return FALSE;
01621 }
01622 /*      END OF FUNCTION atoms_1to4       */
01623 #endif // #ifndef MEM_OPT_VERSION
01624 
01625 //JLai
01626 /************************************************************************/
01627 /*                  */
01628 /*      FUNCTION send_GoMolecule        */
01629 /*                  */
01630 /*  send_Molecule is used by the Master node to distribute the      */
01631 /*   Go information to all the client nodes.  It is NEVER called*/
01632 /*   by the client nodes.              */
01633 /*                  */
01634 /************************************************************************/
01635 void Molecule::send_GoMolecule(MOStream *msg) {
01636   Real *a1, *a2, *a3, *a4;
01637   int *i1, *i2, *i3, *i4;
01638   int maxGoChainsSqr = MAX_GO_CHAINS*MAX_GO_CHAINS;  // JE JLai Go code
01639   msg->put(NumGoChains);
01640   
01641   if (NumGoChains) {
01642     //      int go_indices[MAX_GO_CHAINS+1];        //  Indices from chainIDs to go_array      
01643     //      GoValue go_array[MAX_GO_CHAINS*MAX_GO_CHAINS];   //  Array of Go params
01644     msg->put(MAX_GO_CHAINS+1,go_indices);
01645 
01646     a1 = new Real[maxGoChainsSqr];
01647     a2 = new Real[maxGoChainsSqr];
01648     a3 = new Real[maxGoChainsSqr];
01649     a4 = new Real[maxGoChainsSqr];
01650     i1 = new int[maxGoChainsSqr];
01651     i2 = new int[maxGoChainsSqr];
01652     i3 = new int[maxGoChainsSqr];
01653     i4 = new int[maxGoChainsSqr*MAX_RESTRICTIONS];
01654 
01655     if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) || 
01656          (i1 == NULL) || (i2 == NULL) || (i3 == NULL) || (i4 == NULL) )
01657     {
01658       NAMD_die("memory allocation failed in Molecules::send_Molecules");
01659     }
01660 
01661     for (int i=0; i<maxGoChainsSqr; i++) {
01662       a1[i] = go_array[i].epsilon;
01663       a2[i] = go_array[i].sigmaRep;
01664       a3[i] = go_array[i].epsilonRep;
01665       a4[i] = go_array[i].cutoff;
01666       i1[i] = go_array[i].exp_a;
01667       i2[i] = go_array[i].exp_b;
01668       i3[i] = go_array[i].exp_rep;
01669       for (int j=0; j<MAX_RESTRICTIONS; j++) {
01670         i4[i*MAX_RESTRICTIONS + j] = go_array[i].restrictions[j];
01671       }
01672     }
01673 
01674     msg->put(maxGoChainsSqr, a1);
01675     msg->put(maxGoChainsSqr, a2);
01676     msg->put(maxGoChainsSqr, a3);
01677     msg->put(maxGoChainsSqr, a4);
01678     msg->put(maxGoChainsSqr, i1);
01679     msg->put(maxGoChainsSqr, i2);
01680     msg->put(maxGoChainsSqr, i3);
01681     msg->put(maxGoChainsSqr*MAX_RESTRICTIONS, i4);
01682 
01683     delete [] a1;
01684     delete [] a2;
01685     delete [] a3;
01686     delete [] a4;
01687     delete [] i1;
01688     delete [] i2;
01689     delete [] i3;
01690     delete [] i4;
01691   }
01692 
01693   //Ported JLai
01694   if (simParams->goForcesOn) {
01695     switch(simParams->goMethod) {
01696     case 1:
01697       msg->put(numGoAtoms);
01698       msg->put(numAtoms, goSigmaIndices);
01699       msg->put(numGoAtoms, atomChainTypes);
01700       msg->put(numGoAtoms*numGoAtoms, goSigmas);
01701       msg->put(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
01702       // printf("Molecule.C sending atomChainTypes %d %d \n", numGoAtoms, atomChainTypes);
01703       break;
01704     case 2: //GSS
01705       msg->put(numGoAtoms);
01706       msg->put(numAtoms,pointerToGoBeg);
01707       msg->put(numAtoms,pointerToGoEnd);
01708       msg->put(numAtoms,goSigmaIndices);
01709       msg->put(numAtoms,goResidIndices);
01710       msg->put(numGoAtoms,atomChainTypes);
01711       msg->put(goNumLJPair);
01712       msg->put(goNumLJPair,goIndxLJA);
01713       msg->put(goNumLJPair,goIndxLJB);
01714       msg->put(goNumLJPair,goSigmaPairA);
01715       msg->put(goNumLJPair,goSigmaPairB);
01716       break;
01717     case 3:
01718       msg->put(numGoAtoms);
01719       msg->put(numAtoms, goSigmaIndices);
01720       msg->put(numGoAtoms, atomChainTypes);
01721       //msg->put(numGoAtoms*numGoAtoms, goSigmas);
01722       //msg->put(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
01723       msg->put(numGoAtoms*3, goCoordinates);
01724       msg->put(numGoAtoms, goResids);
01725       break;
01726     }
01727   } 
01728 
01729   msg->end();
01730   delete msg;
01731 }
01732 /*      END OF FUNCTION send_GoMolecule     */
01733 
01734 // JLai
01735 /************************************************************************/
01736 /*                  */
01737 /*      FUNCTION receive_Molecule      */
01738 /*                  */
01739 /*  receive_Molecule is used by all the clients to receive the  */
01740 /*   Go structural data sent out by the master node.  It is NEVER called   */
01741 /*   by the Master node.            */
01742 /*                  */
01743 /************************************************************************/
01744 void Molecule::receive_GoMolecule(MIStream *msg) {
01745       // Ported by JLai -- Original by JE
01746       // JE - receive Go info
01747       Real *a1, *a2, *a3, *a4;
01748       int *i1, *i2, *i3, *i4;
01749       int maxGoChainsSqr = MAX_GO_CHAINS*MAX_GO_CHAINS;  // JE JLai Go code
01750       msg->get(NumGoChains);
01751       
01752       if (NumGoChains) {
01753         //go_indices = new int[MAX_GO_CHAINS+1];
01754         //go_array = new GoValue[MAX_GO_CHAINS*MAX_GO_CHAINS];
01755         
01756         //      int go_indices[MAX_GO_CHAINS+1];        //  Indices from chainIDs to go_array      
01757         //      GoValue go_array[MAX_GO_CHAINS*MAX_GO_CHAINS];   //  Array of Go params
01758         msg->get(MAX_GO_CHAINS+1,go_indices);
01759         
01760         a1 = new Real[maxGoChainsSqr];
01761         a2 = new Real[maxGoChainsSqr];
01762         a3 = new Real[maxGoChainsSqr];
01763         a4 = new Real[maxGoChainsSqr];
01764         i1 = new int[maxGoChainsSqr];
01765         i2 = new int[maxGoChainsSqr];
01766         i3 = new int[maxGoChainsSqr];
01767         i4 = new int[maxGoChainsSqr*MAX_RESTRICTIONS];
01768         
01769         if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) || 
01770              (i1 == NULL) || (i2 == NULL) || (i3 == NULL) || (i4 == NULL) )
01771           {
01772             NAMD_die("memory allocation failed in Molecule::send_Molecule");
01773           }
01774         
01775         msg->get(maxGoChainsSqr, a1);
01776         msg->get(maxGoChainsSqr, a2);
01777         msg->get(maxGoChainsSqr, a3);
01778         msg->get(maxGoChainsSqr, a4);
01779         msg->get(maxGoChainsSqr, i1);
01780         msg->get(maxGoChainsSqr, i2);
01781         msg->get(maxGoChainsSqr, i3);
01782         msg->get(maxGoChainsSqr*MAX_RESTRICTIONS, i4);
01783         
01784         for (int i=0; i<maxGoChainsSqr; i++) {
01785           go_array[i].epsilon = a1[i];
01786           go_array[i].sigmaRep = a2[i];
01787           go_array[i].epsilonRep = a3[i];
01788           go_array[i].cutoff = a4[i];
01789           go_array[i].exp_a = i1[i];
01790           go_array[i].exp_b = i2[i];
01791           go_array[i].exp_rep = i3[i];
01792           for (int j=0; j<MAX_RESTRICTIONS; j++) {
01793             go_array[i].restrictions[j] = i4[i*MAX_RESTRICTIONS + j];
01794           }
01795         }
01796         
01797         delete [] a1;
01798         delete [] a2;
01799         delete [] a3;
01800         delete [] a4;
01801         delete [] i1;
01802         delete [] i2;
01803         delete [] i3;
01804         delete [] i4;
01805 
01806         //msg->get(MAX_GO_CHAINS*MAX_GO_CHAINS, (char*)go_array);
01807         
01808         /*DebugM(3,"NumGoChains:" << NumGoChains << std::endl);
01809           for (int ii=0; ii<MAX_GO_CHAINS; ii++) {
01810           for (int jj=0; jj<MAX_GO_CHAINS; jj++) {
01811           DebugM(3,"go_array[" << ii << "][" << jj << "]:" << go_array[ii][jj] << std::endl);
01812           }
01813           }
01814           for (int ii=0; ii<MAX_GO_CHAINS+1; ii++) {
01815           DebugM(3,"go_indices[" << ii << "]:" << go_indices[ii] << std::endl);
01816           }*/
01817       }
01818 
01819       if (simParams->goForcesOn) {
01820         switch(simParams->goMethod) {
01821         case 1:
01822           msg->get(numGoAtoms);
01823           //printf("Deleting goSigmaIndiciesA\n");
01824           delete [] goSigmaIndices;
01825           goSigmaIndices = new int32[numAtoms];
01826           //printf("Deleting atomChainTypesA\n");
01827           delete [] atomChainTypes;
01828           atomChainTypes = new int32[numGoAtoms];
01829           //printf("Deleting goSigmasA\n");
01830           delete [] goSigmas;
01831           goSigmas = new Real[numGoAtoms*numGoAtoms];
01832           //printf("Deleting goWithinCutoffA\n"); 
01833           delete [] goWithinCutoff;
01834           goWithinCutoff = new bool[numGoAtoms*numGoAtoms];
01835           msg->get(numAtoms, goSigmaIndices);
01836           msg->get(numGoAtoms, atomChainTypes);
01837           msg->get(numGoAtoms*numGoAtoms, goSigmas);
01838           msg->get(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
01839           break;          
01840         case 2: //GSR
01841           msg->get(numGoAtoms);
01842           delete [] pointerToGoBeg;
01843           pointerToGoBeg = new int[numAtoms];
01844           msg->get(numAtoms,pointerToGoBeg);
01845           delete [] pointerToGoEnd;
01846           pointerToGoEnd = new int[numAtoms];
01847           msg->get(numAtoms,pointerToGoEnd);
01848           delete [] goSigmaIndices;
01849           goSigmaIndices = new int32[numAtoms];
01850           msg->get(numAtoms,goSigmaIndices);
01851           delete [] goResidIndices;
01852           goResidIndices = new int32[numAtoms];
01853           msg->get(numAtoms,goResidIndices);      
01854           delete [] atomChainTypes;
01855           atomChainTypes = new int32[numGoAtoms];
01856           msg->get(numGoAtoms,atomChainTypes);
01857           msg->get(goNumLJPair);
01858           delete [] goIndxLJA;
01859           goIndxLJA = new int[goNumLJPair];
01860           msg->get(goNumLJPair,goIndxLJA);
01861           delete [] goIndxLJB;
01862           goIndxLJB = new int[goNumLJPair];
01863           msg->get(goNumLJPair,goIndxLJB);
01864           delete [] goSigmaPairA;
01865           goSigmaPairA = new double[goNumLJPair];
01866           msg->get(goNumLJPair,goSigmaPairA);
01867           delete [] goSigmaPairB;
01868           goSigmaPairB = new double[goNumLJPair];
01869           msg->get(goNumLJPair,goSigmaPairB);  
01870           break;
01871         case 3:
01872           msg->get(numGoAtoms);
01873           //printf("Deleting goSigmaIndiciesB\n");
01874           delete [] goSigmaIndices;
01875           goSigmaIndices = new int32[numAtoms];
01876           //printf("Deleting atomChainTypesB\n");
01877           delete [] atomChainTypes;
01878           atomChainTypes = new int32[numGoAtoms];
01879           //delete [] goSigmas;
01880           //goSigmas = new Real[numGoAtoms*numGoAtoms];
01881           //delete [] goWithinCutoff;
01882           //goWithinCutoff = new bool[numGoAtoms*numGoAtoms];
01883           //printf("Deleting goCoordinatesB\n");
01884           delete [] goCoordinates;
01885           goCoordinates = new Real[numGoAtoms*3];
01886           //printf("Deleting goResidsB\n");
01887           delete [] goResids;
01888           goResids = new int[numGoAtoms];
01889           msg->get(numAtoms, goSigmaIndices);
01890           msg->get(numGoAtoms, atomChainTypes);
01891           //msg->get(numGoAtoms*numGoAtoms, goSigmas);
01892           //msg->get(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
01893           msg->get(numGoAtoms*3, goCoordinates);
01894           msg->get(numGoAtoms, goResids);
01895           break;
01896         }
01897       }
01898 
01899       delete msg;
01900 
01901 }
01902 /*      END OF FUNCTION receive_GoMolecule     */

Generated on Mon Nov 20 01:17:12 2017 for NAMD by  doxygen 1.4.7