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
00026 #include "ComputeGridForce.h"
00027 #include "GridForceGrid.h"
00028 #include "MGridforceParams.h"
00029
00030
00031 #define MIN_DEBUG_LEVEL 3
00032
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
00046
00047
00048
00049
00050
00051
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
00070
00071
00072
00073
00074
00075
00076
00077
00078
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
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113 void Molecule::read_go_file(char *fname)
00114
00115 {
00116
00117 int i;
00118 int j;
00119 int par_type=0;
00120
00121
00122 int skipline;
00123 int skipall = 0;
00124 char buffer[512];
00125 char first_word[512];
00126 int read_count = 0;
00127 int chain1 = 0;
00128 int chain2 = 0;
00129 int int1;
00130 int int2;
00131 Real r1;
00132 char in2[512];
00133 GoValue *goValue1 = NULL;
00134 GoValue *goValue2 = NULL;
00135 Bool sameGoChain = FALSE;
00136 int restrictionCount = 0;
00137 FILE *pfile;
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147 for (i=0; i<MAX_GO_CHAINS+1; i++) {
00148 go_indices[i] = -1;
00149 }
00150
00151
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
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
00175 while (NAMD_read_line(pfile, buffer) != -1)
00176 {
00177
00178 NAMD_find_first_word(buffer, first_word);
00179 skipline=0;
00180
00181
00182
00183
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
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
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;
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
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
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
00327 {
00328
00329
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
00343 fclose(pfile);
00344
00345 return;
00346 }
00347
00348
00349 #if CODE_REDUNDANT
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
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
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
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
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409 Real Molecule::get_go_sigmaRep(int chain1, int chain2)
00410 {
00411 return go_array[MAX_GO_CHAINS*chain1 + chain2].sigmaRep;
00412 }
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
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
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
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
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
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
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
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
00509 #endif
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525 Bool Molecule::go_restricted(int chain1, int chain2, int rDiff)
00526 {
00527 int i;
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
00538
00539
00540
00541
00542
00543
00544
00545
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
00561
00562
00563
00564
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
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;
00582 int bcol = 4;
00583 int32 chainType = 0;
00584 int i;
00585 int j;
00586 int resid1;
00587 int resid2;
00588 int residDiff;
00589 Real sigma;
00590 Real atomAtomDist;
00591 Real exp_a;
00592 Real exp_b;
00593 char filename[NAMD_FILENAME_BUFFER_SIZE];
00594
00595
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
00607
00608
00609 if (goCoordFile == NULL)
00610 {
00611
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
00643 atomChainTypes = new int32[numAtoms];
00644
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
00654 for (i=0; i<numAtoms; i++) {
00655
00656 chainType = (int32)(goPDB->atom(i))->occupancy();
00657
00658 if ( chainType != 0 ) {
00659
00660 atomChainTypes[i] = chainType;
00661 goSigmaIndices[i] = numGoAtoms;
00662 numGoAtoms++;
00663 }
00664 else {
00665 atomChainTypes[i] = 0;
00666 goSigmaIndices[i] = -1;
00667 }
00668
00669 }
00670
00671
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
00681 DebugM(3," numAtoms=" << numAtoms << std::endl);
00682 for (i=0; i<numAtoms; i++) {
00683
00684 resid1 = (goPDB->atom(i))->residueseq();
00685
00686
00687
00688
00689 for (j=i+1; j<numAtoms; j++) {
00690
00691 resid2 = (goPDB->atom(j))->residueseq();
00692
00693
00694
00695
00696
00697
00698 if ( goSigmaIndices[i] != -1 && goSigmaIndices[j] != -1 && !atoms_1to4(i,j) ) {
00699
00700 residDiff = resid2 - resid1;
00701
00702 if (residDiff < 0) residDiff = -residDiff;
00703
00704
00705
00706
00707
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
00739 if (goCoordFile != NULL) {
00740 delete goPDB;
00741 }
00742
00743 return;
00744 }
00745
00746
00747 void Molecule::build_go_sigmas2(StringList *goCoordFile,
00748 char *cwd)
00749 {
00750 DebugM(3,"->build_go_sigmas2" << std::endl);
00751 PDB *goPDB;
00752 int bcol = 4;
00753 int32 chainType = 0;
00754 int32 residType = 0;
00755 int i;
00756 int j;
00757 int resid1;
00758 int resid2;
00759 int residDiff;
00760 Real sigma;
00761 Real atomAtomDist;
00762 Real exp_a;
00763 Real exp_b;
00764 char filename[NAMD_FILENAME_BUFFER_SIZE];
00765
00766
00767 int numLJPair = 0;
00768 long nativeContacts = 0;
00769 long nonnativeContacts = 0;
00770
00771
00772
00773
00774 if (goCoordFile == NULL)
00775 {
00776
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
00808 atomChainTypes = new int32[numAtoms];
00809
00810 goSigmaIndices = new int32[numAtoms];
00811
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
00821 for (i=0; i<numAtoms; i++) {
00822
00823 chainType = (int32)(goPDB->atom(i))->occupancy();
00824
00825 residType = (int32)(goPDB->atom(i))->residueseq();
00826
00827 if ( chainType != 0 ) {
00828
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
00842 ResizeArray<GoPair> tmpGoPair;
00843
00844
00845
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
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
00919 if (goCoordFile != NULL) {
00920 delete goPDB;
00921 }
00922
00923 return;
00924 }
00925
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
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950 void Molecule::build_go_arrays(StringList *goCoordFile,
00951 char *cwd)
00952 {
00953 DebugM(3,"->build_go_arrays" << std::endl);
00954
00955 int bcol = 4;
00956 int32 chainType = 0;
00957 int i;
00958 int j;
00959 BigReal atomAtomDist;
00960 int resid1;
00961 int resid2;
00962 int residDiff;
00963 int goIndex;
00964 int goIndx;
00965 char filename[NAMD_FILENAME_BUFFER_SIZE];
00966
00967
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
00979
00980
00981 if (goCoordFile == NULL)
00982 {
00983
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
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
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
01049
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
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
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
01074
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
01083 energyNative = 0;
01084 energyNonnative = 0;
01085
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
01116 if (goCoordFile != NULL) {
01117 delete goPDB;
01118 }
01119
01120 return;
01121 }
01122
01123 #endif // #ifndef MEM_OPT_VERSION
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134 void Molecule::print_go_sigmas()
01135 {
01136 int i;
01137 int j;
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
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
01162 }
01163 }
01164 if ( goSigmaIndices[i] != -1 ) {
01165 DebugM(3, "-----------" << std::endl);
01166 }
01167 }
01168 return;
01169 }
01170
01171
01172
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
01182 *pairLJEnergy = 0.0;
01183 *pairGaussEnergy = 0.0;
01184
01185
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
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
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
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
01225 }
01226
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
01252 return (LJ + Gauss)*ri;
01253 }
01254 return 0;
01255 }
01256
01257
01258
01259
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
01270
01271 int32 chain1 = atomChainTypes[atom1];
01272 int32 chain2 = atomChainTypes[atom2];
01273
01274
01275 if (chain1 == 0 || chain2 == 0) return 0.0;
01276
01277
01278
01279 Real goCutoff = const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2);
01280
01281 if (goCutoff == 0) return 0.0;
01282
01283
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
01294 }
01295
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
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
01309 *goNative = (4 * epsilon * ( pow1 - pow2 ) );
01310 *goNonnative = 0.0;
01311 }
01312 }
01313 }
01314
01315 return goForce;
01316 }
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
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
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
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
01378
01379
01380
01381
01382
01383
01384
01385
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
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
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
01405
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
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
01414 *goNative = (4 * epsilon * ( pow1 - pow2 ) );
01415 *goNonnative = 0;
01416 }
01417
01418
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
01425 goForce = ((4/(r*r)) * expRep * 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
01428 *goNonnative = (4 * epsilonRep * pow1);
01429 *goNative = 0;
01430 }
01431 }
01432 }
01433
01434
01435 return goForce;
01436 }
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
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
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
01502 return (LJ/r);
01503 } else {
01504
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
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
01521
01522
01523 #ifndef MEM_OPT_VERSION
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537 Bool Molecule::atoms_1to4(unsigned int atom1,
01538 unsigned int atom2)
01539 {
01540 int bondNum;
01541 int angleNum;
01542 int dihedralNum;
01543 int *bonds;
01544 int *angles;
01545 int *dihedrals;
01546 Bond *bond;
01547 Angle *angle;
01548 Dihedral *dihedral;
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
01623 #endif // #ifndef MEM_OPT_VERSION
01624
01625
01626
01627
01628
01629
01630
01631
01632
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;
01639 msg->put(NumGoChains);
01640
01641 if (NumGoChains) {
01642
01643
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
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
01703 break;
01704 case 2:
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
01722
01723 msg->put(numGoAtoms*3, goCoordinates);
01724 msg->put(numGoAtoms, goResids);
01725 break;
01726 }
01727 }
01728
01729 msg->end();
01730 delete msg;
01731 }
01732
01733
01734
01735
01736
01737
01738
01739
01740
01741
01742
01743
01744 void Molecule::receive_GoMolecule(MIStream *msg) {
01745
01746
01747 Real *a1, *a2, *a3, *a4;
01748 int *i1, *i2, *i3, *i4;
01749 int maxGoChainsSqr = MAX_GO_CHAINS*MAX_GO_CHAINS;
01750 msg->get(NumGoChains);
01751
01752 if (NumGoChains) {
01753
01754
01755
01756
01757
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
01807
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817 }
01818
01819 if (simParams->goForcesOn) {
01820 switch(simParams->goMethod) {
01821 case 1:
01822 msg->get(numGoAtoms);
01823
01824 delete [] goSigmaIndices;
01825 goSigmaIndices = new int32[numAtoms];
01826
01827 delete [] atomChainTypes;
01828 atomChainTypes = new int32[numGoAtoms];
01829
01830 delete [] goSigmas;
01831 goSigmas = new Real[numGoAtoms*numGoAtoms];
01832
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:
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
01874 delete [] goSigmaIndices;
01875 goSigmaIndices = new int32[numAtoms];
01876
01877 delete [] atomChainTypes;
01878 atomChainTypes = new int32[numGoAtoms];
01879
01880
01881
01882
01883
01884 delete [] goCoordinates;
01885 goCoordinates = new Real[numGoAtoms*3];
01886
01887 delete [] goResids;
01888 goResids = new int[numGoAtoms];
01889 msg->get(numAtoms, goSigmaIndices);
01890 msg->get(numGoAtoms, atomChainTypes);
01891
01892
01893 msg->get(numGoAtoms*3, goCoordinates);
01894 msg->get(numGoAtoms, goResids);
01895 break;
01896 }
01897 }
01898
01899 delete msg;
01900
01901 }
01902