#include <Molecule.h>
Definition at line 150 of file Molecule.h.
| Molecule::Molecule | ( | SimParameters * | simParams, | |
| Parameters * | param | |||
| ) |
Definition at line 423 of file Molecule.C.
References initialize().
00424 { 00425 initialize(simParams,param); 00426 }
| Molecule::Molecule | ( | SimParameters * | simParams, | |
| Parameters * | param, | |||
| char * | filename, | |||
| ConfigList * | cfgList = NULL | |||
| ) |
Definition at line 436 of file Molecule.C.
References initialize(), SimParameters::LCPOOn, and SimParameters::useCompressedPsf.
00437 { 00438 initialize(simParams,param); 00439 00440 #ifdef MEM_OPT_VERSION 00441 if(simParams->useCompressedPsf) 00442 read_mol_signatures(filename, param, cfgList); 00443 #else 00444 read_psf_file(filename, param); 00445 //LCPO 00446 if (simParams->LCPOOn) 00447 assignLCPOTypes( 0 ); 00448 #endif 00449 }
| Molecule::Molecule | ( | SimParameters * | simParams, | |
| Parameters * | param, | |||
| molfile_plugin_t * | pIOHdl, | |||
| void * | pIOFileHdl, | |||
| int | natoms | |||
| ) |
Definition at line 458 of file Molecule.C.
References initialize(), SimParameters::LCPOOn, NAMD_die(), numAngles, numAtoms, numBonds, numCrossterms, numDihedrals, numImpropers, numRealBonds, setBFactorData(), and setOccupancyData().
00459 { 00460 #ifdef MEM_OPT_VERSION 00461 NAMD_die("Sorry, plugin IO is not supported in the memory optimized version."); 00462 #else 00463 initialize(simParams, param); 00464 numAtoms = natoms; 00465 int optflags = MOLFILE_BADOPTIONS; 00466 molfile_atom_t *atomarray = (molfile_atom_t *) malloc(natoms*sizeof(molfile_atom_t)); 00467 memset(atomarray, 0, natoms*sizeof(molfile_atom_t)); 00468 00469 //1a. read basic atoms information 00470 int rc = pIOHdl->read_structure(pIOFileHdl, &optflags, atomarray); 00471 if (rc != MOLFILE_SUCCESS && rc != MOLFILE_NOSTRUCTUREDATA) { 00472 free(atomarray); 00473 NAMD_die("ERROR: plugin failed reading structure data"); 00474 } 00475 if(optflags == MOLFILE_BADOPTIONS) { 00476 free(atomarray); 00477 NAMD_die("ERROR: plugin didn't initialize optional data flags"); 00478 } 00479 if(optflags & MOLFILE_OCCUPANCY) { 00480 setOccupancyData(atomarray); 00481 } 00482 if(optflags & MOLFILE_BFACTOR) { 00483 setBFactorData(atomarray); 00484 } 00485 //1b. load basic atoms information to the molecule object 00486 plgLoadAtomBasics(atomarray); 00487 free(atomarray); 00488 00489 //2a. read bonds 00490 //indices are one-based in read_bonds 00491 int *from, *to; 00492 float *bondorder; 00493 int *bondtype, nbondtypes; 00494 char **bondtypename; 00495 if(pIOHdl->read_bonds!=NULL) { 00496 if(pIOHdl->read_bonds(pIOFileHdl, &numBonds, &from, &to, &bondorder, 00497 &bondtype, &nbondtypes, &bondtypename)){ 00498 NAMD_die("ERROR: failed reading bond information."); 00499 } 00500 } 00501 //2b. load bonds information to the molecule object 00502 if(numBonds!=0) { 00503 plgLoadBonds(from,to); 00504 } 00505 00506 //3a. read other bonded structures 00507 int *plgAngles, *plgDihedrals, *plgImpropers, *plgCterms; 00508 int ctermcols, ctermrows; 00509 int *angletypes, numangletypes, *dihedraltypes, numdihedraltypes; 00510 int *impropertypes, numimpropertypes; 00511 char **angletypenames, **dihedraltypenames, **impropertypenames; 00512 00513 plgAngles=plgDihedrals=plgImpropers=plgCterms=NULL; 00514 if(pIOHdl->read_angles!=NULL) { 00515 if(pIOHdl->read_angles(pIOFileHdl, 00516 &numAngles, &plgAngles, 00517 &angletypes, &numangletypes, &angletypenames, 00518 &numDihedrals, &plgDihedrals, 00519 &dihedraltypes, &numdihedraltypes, &dihedraltypenames, 00520 &numImpropers, &plgImpropers, 00521 &impropertypes, &numimpropertypes, &impropertypenames, 00522 &numCrossterms, &plgCterms, &ctermcols, &ctermrows)) { 00523 NAMD_die("ERROR: failed reading angle information."); 00524 } 00525 } 00526 //3b. load other bonded structures to the molecule object 00527 if(numAngles!=0) plgLoadAngles(plgAngles); 00528 if(numDihedrals!=0) plgLoadDihedrals(plgDihedrals); 00529 if(numImpropers!=0) plgLoadImpropers(plgImpropers); 00530 if(numCrossterms!=0) plgLoadCrossterms(plgCterms); 00531 00532 numRealBonds = numBonds; 00533 build_atom_status(); 00534 //LCPO 00535 if (simParams->LCPOOn) 00536 assignLCPOTypes( 2 ); 00537 #endif 00538 }
| Molecule::Molecule | ( | SimParameters * | , | |
| Parameters * | , | |||
| Ambertoppar * | ||||
| ) |
| Molecule::Molecule | ( | SimParameters * | , | |
| Parameters * | , | |||
| const GromacsTopFile * | ||||
| ) |
| Molecule::~Molecule | ( | ) |
Definition at line 552 of file Molecule.C.
References atomSigPool, numAtoms, ss_index, and ss_vdw_type.
00553 { 00554 /* Check to see if each array was ever allocated. If it was */ 00555 /* then free it */ 00556 if (atoms != NULL) 00557 delete [] atoms; 00558 00559 if (atomNames != NULL) 00560 { 00561 // subarrarys allocated from arena - automatically deleted 00562 delete [] atomNames; 00563 } 00564 delete nameArena; 00565 00566 if (resLookup != NULL) 00567 delete resLookup; 00568 00569 // DRUDE: free arrays read from PSF 00570 if (drudeConsts != NULL) delete [] drudeConsts; 00571 if (lphosts != NULL) delete [] lphosts; 00572 if (anisos != NULL) delete [] anisos; 00573 if (tholes != NULL) delete [] tholes; 00574 if (lphostIndexes != NULL) delete [] lphostIndexes; 00575 // DRUDE 00576 00577 //LCPO 00578 if (lcpoParamType != NULL) delete [] lcpoParamType; 00579 00580 #ifdef MEM_OPT_VERSION 00581 if(eachAtomSig) delete [] eachAtomSig; 00582 if(atomSigPool) delete [] atomSigPool; 00583 #else 00584 if (bonds != NULL) 00585 delete [] bonds; 00586 00587 if (angles != NULL) 00588 delete [] angles; 00589 00590 if (dihedrals != NULL) 00591 delete [] dihedrals; 00592 00593 if (impropers != NULL) 00594 delete [] impropers; 00595 00596 if (crossterms != NULL) 00597 delete [] crossterms; 00598 00599 if (exclusions != NULL) 00600 delete [] exclusions; 00601 #endif 00602 00603 if (donors != NULL) 00604 delete [] donors; 00605 00606 if (acceptors != NULL) 00607 delete [] acceptors; 00608 00609 #ifdef MEM_OPT_VERSION 00610 if(exclSigPool) delete [] exclSigPool; 00611 if(exclChkSigPool) delete [] exclChkSigPool; 00612 if(eachAtomExclSig) delete [] eachAtomExclSig; 00613 if(fixedAtomsSet) delete fixedAtomsSet; 00614 if(constrainedAtomsSet) delete constrainedAtomsSet; 00615 #else 00616 if (bondsByAtom != NULL) 00617 delete [] bondsByAtom; 00618 00619 if (anglesByAtom != NULL) 00620 delete [] anglesByAtom; 00621 00622 if (dihedralsByAtom != NULL) 00623 delete [] dihedralsByAtom; 00624 00625 if (impropersByAtom != NULL) 00626 delete [] impropersByAtom; 00627 00628 if (crosstermsByAtom != NULL) 00629 delete [] crosstermsByAtom; 00630 00631 if (exclusionsByAtom != NULL) 00632 delete [] exclusionsByAtom; 00633 00634 if (fullExclusionsByAtom != NULL) 00635 delete [] fullExclusionsByAtom; 00636 00637 if (modExclusionsByAtom != NULL) 00638 delete [] modExclusionsByAtom; 00639 00640 if (all_exclusions != NULL) 00641 delete [] all_exclusions; 00642 00643 // JLai 00644 if (gromacsPairByAtom != NULL) 00645 delete [] gromacsPairByAtom; 00646 // End of JLai 00647 00648 // DRUDE 00649 if (tholesByAtom != NULL) 00650 delete [] tholesByAtom; 00651 if (anisosByAtom != NULL) 00652 delete [] anisosByAtom; 00653 // DRUDE 00654 #endif 00655 00656 //LCPO 00657 if (lcpoParamType != NULL) 00658 delete [] lcpoParamType; 00659 00660 if (fixedAtomFlags != NULL) 00661 delete [] fixedAtomFlags; 00662 00663 if (stirIndexes != NULL) 00664 delete [] stirIndexes; 00665 00666 00667 #ifdef MEM_OPT_VERSION 00668 if(clusterSigs != NULL){ 00669 delete [] clusterSigs; 00670 } 00671 #else 00672 if (cluster != NULL) 00673 delete [] cluster; 00674 #endif 00675 if (clusterSize != NULL) 00676 delete [] clusterSize; 00677 00678 if (exPressureAtomFlags != NULL) 00679 delete [] exPressureAtomFlags; 00680 00681 if (rigidBondLengths != NULL) 00682 delete [] rigidBondLengths; 00683 00684 //fepb 00685 if (fepAtomFlags != NULL) 00686 delete [] fepAtomFlags; 00687 //fepe 00688 00689 //soluteScaling 00690 if (ssAtomFlags != NULL) 00691 delete [] ssAtomFlags; 00692 if (ss_vdw_type != NULL) 00693 delete [] ss_vdw_type; 00694 if (ss_index != NULL) 00695 delete [] ss_index; 00696 //soluteScaling 00697 00698 if (qmAtomGroup != NULL) 00699 delete [] qmAtomGroup; 00700 00701 if (qmAtmIndx != NULL) 00702 delete [] qmAtmIndx; 00703 00704 if (qmAtmChrg != NULL) 00705 delete [] qmAtmChrg; 00706 00707 00708 if (qmGrpNumBonds != NULL) 00709 delete [] qmGrpNumBonds; 00710 00711 if (qmGrpSizes != NULL) 00712 delete [] qmGrpSizes; 00713 00714 if (qmDummyBondVal != NULL) 00715 delete [] qmDummyBondVal; 00716 00717 if (qmMMNumTargs != NULL) 00718 delete [] qmMMNumTargs; 00719 00720 if (qmGrpID != NULL) 00721 delete [] qmGrpID; 00722 00723 if (qmGrpChrg != NULL) 00724 delete [] qmGrpChrg; 00725 00726 if (qmGrpMult != NULL) 00727 delete [] qmGrpMult; 00728 00729 if (qmMeMMindx != NULL) 00730 delete [] qmMeMMindx; 00731 00732 if (qmMeQMGrp != NULL) 00733 delete [] qmMeQMGrp; 00734 00735 if (qmLSSSize != NULL) 00736 delete [] qmLSSSize; 00737 00738 if (qmLSSIdxs != NULL) 00739 delete [] qmLSSIdxs; 00740 00741 if (qmLSSMass != NULL) 00742 delete [] qmLSSMass; 00743 00744 if (qmLSSRefSize != NULL) 00745 delete [] qmLSSRefSize; 00746 00747 if (qmLSSRefIDs != NULL) 00748 delete [] qmLSSRefIDs; 00749 00750 if (qmLSSRefMass != NULL) 00751 delete [] qmLSSRefMass; 00752 00753 if (qmMMBond != NULL) { 00754 for(int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) { 00755 if (qmMMBond[grpIndx] != NULL) 00756 delete [] qmMMBond[grpIndx]; 00757 } 00758 delete [] qmMMBond; 00759 } 00760 00761 if (qmGrpBonds != NULL) { 00762 for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) { 00763 if (qmGrpBonds[grpIndx] != NULL) 00764 delete [] qmGrpBonds[grpIndx]; 00765 } 00766 delete [] qmGrpBonds; 00767 } 00768 00769 if (qmMMBondedIndx != NULL) { 00770 for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) { 00771 if (qmMMBondedIndx[grpIndx] != NULL) 00772 delete [] qmMMBondedIndx[grpIndx]; 00773 } 00774 delete [] qmMMBondedIndx; 00775 } 00776 00777 if (qmMMChargeTarget != NULL) { 00778 for(int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) { 00779 if (qmMMChargeTarget[grpIndx] != NULL) 00780 delete [] qmMMChargeTarget[grpIndx]; 00781 } 00782 delete [] qmMMChargeTarget; 00783 } 00784 00785 if (qmElementArray != NULL){ 00786 for(int atmInd = 0 ; atmInd < numAtoms; atmInd++) { 00787 if (qmElementArray[atmInd] != NULL) 00788 delete [] qmElementArray[atmInd]; 00789 } 00790 delete [] qmElementArray; 00791 } 00792 00793 if (qmDummyElement != NULL){ 00794 for(int atmInd = 0 ; atmInd < numAtoms; atmInd++) { 00795 if (qmDummyElement[atmInd] != NULL) 00796 delete [] qmDummyElement[atmInd]; 00797 } 00798 delete [] qmDummyElement; 00799 } 00800 00801 if (qmCustPCSizes != NULL){ 00802 delete [] qmCustPCSizes; 00803 } 00804 00805 if (qmCustomPCIdxs != NULL){ 00806 delete [] qmCustomPCIdxs; 00807 } 00808 00809 if (cSMDindex != NULL) { 00810 for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) { 00811 if (cSMDindex[grpIndx] != NULL) 00812 delete [] cSMDindex[grpIndx]; 00813 } 00814 delete [] cSMDindex; 00815 } 00816 00817 if (cSMDpairs != NULL) { 00818 for(int instIndx = 0 ; instIndx < cSMDnumInst; instIndx++) { 00819 if (cSMDpairs[instIndx] != NULL) 00820 delete [] cSMDpairs[instIndx]; 00821 } 00822 delete [] cSMDpairs; 00823 } 00824 00825 if (cSMDindxLen != NULL) 00826 delete [] cSMDindxLen; 00827 if (cSMDKs != NULL) 00828 delete [] cSMDKs; 00829 if (cSMDVels != NULL) 00830 delete [] cSMDVels; 00831 if (cSMDcoffs != NULL) 00832 delete [] cSMDcoffs; 00833 00834 #ifndef MEM_OPT_VERSION 00835 delete arena; 00836 delete exclArena; 00837 #endif 00838 }
| Real Molecule::atomcharge | ( | int | anum | ) | const [inline] |
Definition at line 1035 of file Molecule.h.
References atoms.
Referenced by WorkDistrib::createAtomLists(), NamdState::loadStructure(), Sequencer::reloadCharges(), Sequencer::rescaleSoluteCharges(), colvarproxy_namd::update_atom_properties(), and colvarproxy_namd::update_group_properties().
| Real Molecule::atommass | ( | int | anum | ) | const [inline] |
Definition at line 1025 of file Molecule.h.
References atoms.
Referenced by WorkDistrib::createAtomLists(), GlobalMasterFreeEnergy::getMass(), GlobalMasterEasy::getMass(), NamdState::loadStructure(), Tcl_centerOfMass(), Tcl_radiusOfGyration(), colvarproxy_namd::update_atom_properties(), and colvarproxy_namd::update_group_properties().
| Bool Molecule::atoms_1to4 | ( | unsigned int | atom1, | |
| unsigned int | atom2 | |||
| ) |
Definition at line 1537 of file GoMolecule.C.
References dihedral::atom1, angle::atom1, bond::atom1, dihedral::atom2, angle::atom2, bond::atom2, dihedral::atom3, angle::atom3, dihedral::atom4, DebugM, FALSE, get_angle(), get_angles_for_atom(), get_bond(), get_bonds_for_atom(), get_dihedral(), get_dihedrals_for_atom(), and TRUE.
Referenced by build_go_arrays(), build_go_sigmas(), and build_go_sigmas2().
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 }
| Index Molecule::atomvdwtype | ( | int | anum | ) | const [inline] |
Definition at line 1045 of file Molecule.h.
References atoms.
Referenced by WorkDistrib::createAtomLists(), and dumpbench().
| void Molecule::build_constant_forces | ( | char * | ) |
Referenced by NamdState::loadStructure().
| void Molecule::build_constorque_params | ( | StringList * | , | |
| StringList * | , | |||
| StringList * | , | |||
| StringList * | , | |||
| StringList * | , | |||
| StringList * | , | |||
| PDB * | , | |||
| char * | ||||
| ) |
Referenced by NamdState::loadStructure().
| void Molecule::build_constraint_params | ( | StringList * | , | |
| StringList * | , | |||
| StringList * | , | |||
| PDB * | , | |||
| char * | ||||
| ) |
Referenced by NamdState::loadStructure().
| void Molecule::build_exPressure_atoms | ( | StringList * | , | |
| StringList * | , | |||
| PDB * | , | |||
| char * | ||||
| ) |
Referenced by NamdState::loadStructure().
| void Molecule::build_extra_bonds | ( | Parameters * | parameters, | |
| StringList * | file | |||
| ) |
Referenced by NamdState::loadStructure().
| void Molecule::build_fep_flags | ( | StringList * | , | |
| StringList * | , | |||
| PDB * | , | |||
| char * | , | |||
| const char * | ||||
| ) |
Referenced by NamdState::loadStructure().
| void Molecule::build_fixed_atoms | ( | StringList * | , | |
| StringList * | , | |||
| PDB * | , | |||
| char * | ||||
| ) |
Referenced by NamdState::loadStructure().
| void Molecule::build_go_arrays | ( | StringList * | goCoordFile, | |
| char * | cwd | |||
| ) |
goSigmas = new Real[numGoAtoms*numGoAtoms]; goWithinCutoff = new bool[numGoAtoms*numGoAtoms]; for (i=0; i<numGoAtoms; i++) { for (j=0; j<numGoAtoms; j++) { goSigmas[i*numGoAtoms + j] = 0.0; goWithinCutoff[i*numGoAtoms + j] = false; } }
Definition at line 950 of file GoMolecule.C.
References PDB::atom(), atomChainTypes, atoms_1to4(), StringList::data, DebugM, endi(), energyNative, energyNonnative, get_go_cutoff(), go_restricted(), goCoordinates, goPDB, goResids, goSigmaIndices, iINFO(), iout, j, NAMD_die(), StringList::next, PDB::num_atoms(), numAtoms, numGoAtoms, PDBAtom::residueseq(), PDBAtom::xcoor(), PDBAtom::ycoor(), and PDBAtom::zcoor().
Referenced by NamdState::loadStructure().
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 }
| void Molecule::build_go_params | ( | StringList * | g | ) |
Definition at line 79 of file GoMolecule.C.
References StringList::data, endi(), iINFO(), iout, NAMD_die(), StringList::next, and read_go_file().
Referenced by NamdState::loadStructure().
00079 { 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 }
| void Molecule::build_go_sigmas | ( | StringList * | goCoordFile, | |
| char * | cwd | |||
| ) |
Definition at line 577 of file GoMolecule.C.
References PDB::atom(), atomChainTypes, atoms_1to4(), StringList::data, DebugM, endi(), get_go_cutoff(), get_go_exp_a(), get_go_exp_b(), go_restricted(), goPDB, goSigmaIndices, goSigmas, goWithinCutoff, iINFO(), iout, j, NAMD_die(), StringList::next, PDB::num_atoms(), numAtoms, and numGoAtoms.
Referenced by NamdState::loadStructure().
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 }
| void Molecule::build_go_sigmas2 | ( | StringList * | goCoordFile, | |
| char * | cwd | |||
| ) |
Definition at line 747 of file GoMolecule.C.
References go_pair::A, ResizeArray< Elem >::add(), PDB::atom(), atomChainTypes, atoms_1to4(), go_pair::B, ResizeArray< Elem >::begin(), StringList::data, DebugM, ResizeArray< Elem >::end(), endi(), get_go_cutoff(), get_go_exp_a(), get_go_exp_b(), go_restricted(), go_pair::goIndxA, go_pair::goIndxB, goIndxLJA, goIndxLJB, goNumLJPair, goPDB, goResidIndices, goSigmaIndices, goSigmaPairA, goSigmaPairB, iINFO(), iout, j, NAMD_die(), StringList::next, PDB::num_atoms(), numAtoms, numGoAtoms, numLJPair, pointerToGoBeg, pointerToGoEnd, and sort.
Referenced by NamdState::loadStructure().
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 }
| void Molecule::build_gridforce_params | ( | StringList * | gridfrcfile, | |
| StringList * | gridfrccol, | |||
| StringList * | gridfrcchrgcol, | |||
| StringList * | potfile, | |||
| PDB * | initial_pdb, | |||
| char * | cwd | |||
| ) |
Definition at line 6171 of file Molecule.C.
References PDB::atom(), DebugM, endi(), MGridforceParamsList::get_first(), GridforceGrid::get_total_grids(), MGridforceParams::gridforceCol, MGridforceParams::gridforceFile, MGridforceParams::gridforceQcol, MGridforceParams::gridforceVfile, iout, iWARN(), SimParameters::mgridforcelist, NAMD_die(), GridforceGrid::new_grid(), MGridforceParams::next, PDB::num_atoms(), numGridforceGrids, and numGridforces.
Referenced by NamdState::loadStructure().
06177 { 06178 PDB *kPDB; 06179 register int i; // Loop counters 06180 register int j; 06181 register int k; 06182 06183 DebugM(3, "Entered build_gridforce_params multi...\n"); 06184 // DebugM(3, "\tgridfrcfile = " << gridfrcfile->data << endi); 06185 // DebugM(3, "\tgridfrccol = " << gridfrccol->data << endi); 06186 06187 MGridforceParams* mgridParams = simParams->mgridforcelist.get_first(); 06188 numGridforceGrids = 0; 06189 while (mgridParams != NULL) { 06190 numGridforceGrids++; 06191 mgridParams = mgridParams->next; 06192 } 06193 06194 DebugM(3, "numGridforceGrids = " << numGridforceGrids << "\n"); 06195 gridfrcIndexes = new int32*[numGridforceGrids]; 06196 gridfrcParams = new GridforceParams*[numGridforceGrids]; 06197 gridfrcGrid = new GridforceGrid*[numGridforceGrids]; 06198 numGridforces = new int[numGridforceGrids]; 06199 06200 int grandTotalGrids = 0; // including all subgrids 06201 06202 mgridParams = simParams->mgridforcelist.get_first(); 06203 for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) { 06204 int current_index=0; // Index into values used 06205 int kcol = 5; // Column to look for force constant in 06206 int qcol = 0; // Column for charge (default 0: use electric charge) 06207 Real kval = 0; // Force constant value retreived 06208 char filename[129]; // PDB filename 06209 char potfilename[129]; // Potential file name 06210 06211 if (mgridParams == NULL) { 06212 NAMD_die("Problem with mgridParams!"); 06213 } 06214 06215 // Now load values from mgridforcelist object 06216 if (mgridParams->gridforceFile == NULL) 06217 { 06218 if ( ! initial_pdb ) NAMD_die("Initial PDB file unavailable, gridforceFile required."); 06219 kPDB = initial_pdb; 06220 } 06221 else 06222 { 06223 DebugM(4, "mgridParams->gridforceFile = " << mgridParams->gridforceFile << "\n" << endi); 06224 06225 if ( (cwd == NULL) || (mgridParams->gridforceFile[0] == '/') ) 06226 { 06227 strcpy(filename, mgridParams->gridforceFile); 06228 } 06229 else 06230 { 06231 strcpy(filename, cwd); 06232 strcat(filename, mgridParams->gridforceFile); 06233 } 06234 06235 kPDB = new PDB(filename); 06236 if ( kPDB == NULL ) 06237 { 06238 NAMD_die("Memory allocation failed in Molecule::build_gridforce_params"); 06239 } 06240 06241 if (kPDB->num_atoms() != numAtoms) 06242 { 06243 NAMD_die("Number of atoms in grid force PDB doesn't match coordinate PDB"); 06244 } 06245 } 06246 06247 // Get the column that the force constant is going to be in. It 06248 // can be in any of the 5 floating point fields in the PDB, according 06249 // to what the user wants. The allowable fields are X, Y, Z, O, or 06250 // B which correspond to the 1st, 2nd, ... 5th floating point fields. 06251 // The default is the 5th field, which is beta (temperature factor) 06252 if (mgridParams->gridforceCol == NULL) 06253 { 06254 kcol = 5; 06255 } 06256 else 06257 { 06258 if (strcasecmp(mgridParams->gridforceCol, "X") == 0) 06259 { 06260 kcol=1; 06261 } 06262 else if (strcasecmp(mgridParams->gridforceCol, "Y") == 0) 06263 { 06264 kcol=2; 06265 } 06266 else if (strcasecmp(mgridParams->gridforceCol, "Z") == 0) 06267 { 06268 kcol=3; 06269 } 06270 else if (strcasecmp(mgridParams->gridforceCol, "O") == 0) 06271 { 06272 kcol=4; 06273 } 06274 else if (strcasecmp(mgridParams->gridforceCol, "B") == 0) 06275 { 06276 kcol=5; 06277 } 06278 else 06279 { 06280 NAMD_die("gridforcecol must have value of X, Y, Z, O, or B"); 06281 } 06282 } 06283 06284 // Get the column that the charge is going to be in. 06285 if (mgridParams->gridforceQcol == NULL) 06286 { 06287 qcol = 0; // Default: don't read charge from file, use electric charge 06288 } 06289 else 06290 { 06291 if (strcasecmp(mgridParams->gridforceQcol, "X") == 0) 06292 { 06293 qcol=1; 06294 } 06295 else if (strcasecmp(mgridParams->gridforceQcol, "Y") == 0) 06296 { 06297 qcol=2; 06298 } 06299 else if (strcasecmp(mgridParams->gridforceQcol, "Z") == 0) 06300 { 06301 qcol=3; 06302 } 06303 else if (strcasecmp(mgridParams->gridforceQcol, "O") == 0) 06304 { 06305 qcol=4; 06306 } 06307 else if (strcasecmp(mgridParams->gridforceQcol, "B") == 0) 06308 { 06309 qcol=5; 06310 } 06311 else 06312 { 06313 NAMD_die("gridforcechargecol must have value of X, Y, Z, O, or B"); 06314 } 06315 } 06316 06317 if (kcol == qcol) { 06318 NAMD_die("gridforcecol and gridforcechargecol cannot have same value"); 06319 } 06320 06321 06322 // Allocate an array that will store an index into the constraint 06323 // parameters for each atom. If the atom is not constrained, its 06324 // value will be set to -1 in this array. 06325 gridfrcIndexes[gridnum] = new int32[numAtoms]; 06326 06327 if (gridfrcIndexes[gridnum] == NULL) 06328 { 06329 NAMD_die("memory allocation failed in Molecule::build_gridforce_params()"); 06330 } 06331 06332 // Loop through all the atoms and find out which ones are constrained 06333 for (i=0; i<numAtoms; i++) 06334 { 06335 // Get the k value based on where we were told to find it 06336 switch (kcol) 06337 { 06338 case 1: 06339 kval = (kPDB->atom(i))->xcoor(); 06340 break; 06341 case 2: 06342 kval = (kPDB->atom(i))->ycoor(); 06343 break; 06344 case 3: 06345 kval = (kPDB->atom(i))->zcoor(); 06346 break; 06347 case 4: 06348 kval = (kPDB->atom(i))->occupancy(); 06349 break; 06350 case 5: 06351 kval = (kPDB->atom(i))->temperaturefactor(); 06352 break; 06353 } 06354 06355 if (kval > 0.0) 06356 { 06357 // This atom is constrained 06358 gridfrcIndexes[gridnum][i] = current_index; 06359 current_index++; 06360 } 06361 else 06362 { 06363 // This atom is not constrained 06364 gridfrcIndexes[gridnum][i] = -1; 06365 } 06366 } 06367 06368 if (current_index == 0) 06369 { 06370 // Constraints were turned on, but there weren't really any constrained 06371 iout << iWARN << "NO GRIDFORCE ATOMS WERE FOUND, BUT GRIDFORCE IS ON . . .\n" << endi; 06372 } 06373 else 06374 { 06375 // Allocate an array to hold the constraint parameters 06376 gridfrcParams[gridnum] = new GridforceParams[current_index]; 06377 if (gridfrcParams[gridnum] == NULL) 06378 { 06379 NAMD_die("memory allocation failed in Molecule::build_gridforce_params"); 06380 } 06381 } 06382 06383 numGridforces[gridnum] = current_index; 06384 06385 // Loop through all the atoms and assign the parameters for those 06386 // that are constrained 06387 for (i=0; i<numAtoms; i++) 06388 { 06389 if (gridfrcIndexes[gridnum][i] != -1) 06390 { 06391 // This atom has grid force, so get the k value again 06392 switch (kcol) 06393 { 06394 case 1: 06395 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->xcoor(); 06396 break; 06397 case 2: 06398 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->ycoor(); 06399 break; 06400 case 3: 06401 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->zcoor(); 06402 break; 06403 case 4: 06404 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->occupancy(); 06405 break; 06406 case 5: 06407 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->temperaturefactor(); 06408 break; 06409 } 06410 06411 // Also get charge column 06412 switch (qcol) 06413 { 06414 case 0: 06415 #ifdef MEM_OPT_VERSION 06416 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atomChargePool[eachAtomCharge[i]]; 06417 #else 06418 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atoms[i].charge; 06419 #endif 06420 break; 06421 case 1: 06422 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->xcoor(); 06423 break; 06424 case 2: 06425 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->ycoor(); 06426 break; 06427 case 3: 06428 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->zcoor(); 06429 break; 06430 case 4: 06431 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->occupancy(); 06432 break; 06433 case 5: 06434 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->temperaturefactor(); 06435 break; 06436 } 06437 } 06438 } 06439 06440 // If we had to create new PDB objects, delete them now 06441 if (mgridParams->gridforceFile != NULL) 06442 { 06443 delete kPDB; 06444 } 06445 06446 // Now we fill in our grid information 06447 06448 // Open potential file 06449 if ( (cwd == NULL) || (mgridParams->gridforceVfile[0] == '/') ) 06450 { 06451 strcpy(potfilename, mgridParams->gridforceVfile); 06452 } 06453 else 06454 { 06455 strcpy(potfilename, cwd); 06456 strcat(potfilename, mgridParams->gridforceVfile); 06457 } 06458 06459 // iout << iINFO << "Allocating grid " << gridnum 06460 // << "\n" << endi; 06461 06462 DebugM(3, "allocating GridforceGrid(" << gridnum << ")\n"); 06463 gridfrcGrid[gridnum] = GridforceGrid::new_grid(gridnum, potfilename, simParams, mgridParams); 06464 06465 grandTotalGrids += gridfrcGrid[gridnum]->get_total_grids(); 06466 DebugM(4, "grandTotalGrids = " << grandTotalGrids << "\n" << endi); 06467 06468 // Finally, get next mgridParams pointer 06469 mgridParams = mgridParams->next; 06470 } 06471 }
| void Molecule::build_gro_pair | ( | ) |
| void Molecule::build_langevin_params | ( | StringList * | , | |
| StringList * | , | |||
| PDB * | , | |||
| char * | ||||
| ) |
Referenced by NamdState::loadStructure().
| void Molecule::build_movdrag_params | ( | StringList * | , | |
| StringList * | , | |||
| StringList * | , | |||
| PDB * | , | |||
| char * | ||||
| ) |
Referenced by NamdState::loadStructure().
| void Molecule::build_rotdrag_params | ( | StringList * | , | |
| StringList * | , | |||
| StringList * | , | |||
| StringList * | , | |||
| StringList * | , | |||
| StringList * | , | |||
| PDB * | , | |||
| char * | ||||
| ) |
Referenced by NamdState::loadStructure().
| void Molecule::build_ss_flags | ( | const StringList * | ssfile, | |
| const StringList * | sscol, | |||
| PDB * | initial_pdb, | |||
| const char * | cwd | |||
| ) |
Build the flags needed for solute scaling.
A PDB file is read, indicating which atoms are to be scaled, and an array is maintained marking which are to be included. Each marked atom then has its corresponding van der Waals type number reassigned to enable extending the LJTable with scaled interaction values.
| ssfile | config "soluteScalingFile = my.pdb" for PDB filename | |
| sscol | config "soluteScalingCol = O" for column of PDB ATOM records | |
| initial_pdb | the initial PDB file | |
| cwd | current working directory |
Referenced by NamdState::loadStructure().
| void Molecule::build_stirred_atoms | ( | StringList * | , | |
| StringList * | , | |||
| PDB * | , | |||
| char * | ||||
| ) |
Referenced by NamdState::loadStructure().
| int Molecule::checkexcl | ( | int | atom1, | |
| int | atom2 | |||
| ) | const |
| void Molecule::compute_LJcorrection | ( | ) |
Referenced by NamdState::loadStructure().
| void Molecule::delete_alch_bonded | ( | void | ) |
Referenced by NamdState::loadStructure().
| void Molecule::delete_qm_bonded | ( | void | ) |
Definition at line 1547 of file MoleculeQM.C.
References crossterm::atom1, improper::atom1, dihedral::atom1, angle::atom1, bond::atom1, DebugM, endi(), SimParameters::extraBondsOn, iERROR(), iINFO(), iout, NAMD_die(), numAngles, numBonds, numCrossterms, numDihedrals, numImpropers, numRealBonds, and ResizeArray< Elem >::size().
01547 { 01548 01549 #ifdef MEM_OPT_VERSION 01550 NAMD_die("QMMM interface is not supported in memory optimized builds."); 01551 #else 01552 01553 DebugM(3,"Cleaning QM bonds, angles, dihedrals, impropers and crossterms.\n") 01554 01555 // Bonds 01556 Bond *nonQMBonds; 01557 nonQMBonds = new Bond[numBonds] ; 01558 int nonQMBondCount = 0; 01559 qmDroppedBonds = 0; 01560 for (int i = 0; i < numBonds; i++) { 01561 01562 int part1 = qmAtomGroup[bonds[i].atom1]; 01563 int part2 = qmAtomGroup[bonds[i].atom2]; 01564 if (part1 > 0 && part2 > 0 ) { 01565 01566 // If the user defined extra bonds, we check if the QM-QM bond is an "extra" 01567 // bond, and do not delete it. 01568 if (simParams->extraBondsOn) { 01569 // std::cout << "Checking Bond: " << bonds[i].atom1 << "," << bonds[i].atom2 << std::endl ; 01570 01571 for (int ebi=0; ebi < qmExtraBonds.size() ; ebi++) { 01572 01573 if( (qmExtraBonds[ebi].atom1 == bonds[i].atom1 || 01574 qmExtraBonds[ebi].atom1 == bonds[i].atom2) && 01575 (qmExtraBonds[ebi].atom2 == bonds[i].atom1 || 01576 qmExtraBonds[ebi].atom2 == bonds[i].atom2) ) { 01577 // std::cout << "This is an extra bond! We will keep it." << std::endl; 01578 nonQMBonds[nonQMBondCount++] = bonds[i]; 01579 break; 01580 } 01581 } 01582 } 01583 01584 qmDroppedBonds++; 01585 } else { 01586 // Just a simple sanity check. 01587 // If there are no QM bonds defined by the user but we find a 01588 // bond between a QM and an MM atom, error out. 01589 if ((part1 > 0 || part2 > 0) && qmNumBonds == 0) { 01590 iout << iERROR << " Atoms " << bonds[i].atom1 01591 << " and " << bonds[i].atom2 << " form a QM-MM bond!\n" << endi ; 01592 NAMD_die("No QM-MM bonds were defined by the user, but one was found."); 01593 } 01594 nonQMBonds[nonQMBondCount++] = bonds[i]; 01595 } 01596 } 01597 numBonds = nonQMBondCount; 01598 delete [] bonds; 01599 bonds = new Bond[numBonds]; 01600 for (int i = 0; i < nonQMBondCount; i++) { 01601 bonds[i]=nonQMBonds[i]; 01602 } 01603 delete [] nonQMBonds; 01604 numRealBonds = numBonds; 01605 01606 // Angles 01607 Angle *nonQMAngles; 01608 nonQMAngles = new Angle[numAngles]; 01609 int nonQMAngleCount = 0; 01610 qmDroppedAngles = 0; 01611 for (int i = 0; i < numAngles; i++) { 01612 int part1 = qmAtomGroup[angles[i].atom1]; 01613 int part2 = qmAtomGroup[angles[i].atom2]; 01614 int part3 = qmAtomGroup[angles[i].atom3]; 01615 if (part1 > 0 && part2 > 0 && part3 > 0) { 01616 //CkPrintf("-----ANGLE ATOMS %i %i %i partitions %i %i %i\n",angles[i].atom1, angles[i].atom2, angles[i].atom3, part1, part2, part3); 01617 qmDroppedAngles++; 01618 } 01619 else { 01620 nonQMAngles[nonQMAngleCount++] = angles[i]; 01621 } 01622 } 01623 numAngles = nonQMAngleCount; 01624 delete [] angles; 01625 angles = new Angle[numAngles]; 01626 for (int i = 0; i < nonQMAngleCount; i++) { 01627 angles[i]=nonQMAngles[i]; 01628 } 01629 delete [] nonQMAngles; 01630 01631 01632 // Dihedrals 01633 Dihedral *nonQMDihedrals; 01634 nonQMDihedrals = new Dihedral[numDihedrals]; 01635 int nonQMDihedralCount = 0; 01636 qmDroppedDihedrals = 0; 01637 for (int i = 0; i < numDihedrals; i++) { 01638 int part1 = qmAtomGroup[dihedrals[i].atom1]; 01639 int part2 = qmAtomGroup[dihedrals[i].atom2]; 01640 int part3 = qmAtomGroup[dihedrals[i].atom3]; 01641 int part4 = qmAtomGroup[dihedrals[i].atom4]; 01642 if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0) { 01643 //CkPrintf("-----i %i DIHEDRAL ATOMS %i %i %i %i partitions %i %i %i %i\n",i,dihedrals[i].atom1, dihedrals[i].atom2, dihedrals[i].atom3, dihedrals[i].atom4, part1, part2, part3,part4); 01644 qmDroppedDihedrals++; 01645 } 01646 else { 01647 nonQMDihedrals[nonQMDihedralCount++] = dihedrals[i]; 01648 } 01649 } 01650 numDihedrals = nonQMDihedralCount; 01651 delete [] dihedrals; 01652 dihedrals = new Dihedral[numDihedrals]; 01653 for (int i = 0; i < numDihedrals; i++) { 01654 dihedrals[i]=nonQMDihedrals[i]; 01655 } 01656 delete [] nonQMDihedrals; 01657 01658 // Impropers 01659 Improper *nonQMImpropers; 01660 nonQMImpropers = new Improper[numImpropers]; 01661 int nonQMImproperCount = 0; 01662 qmDroppedImpropers = 0; 01663 for (int i = 0; i < numImpropers; i++) { 01664 int part1 = qmAtomGroup[impropers[i].atom1]; 01665 int part2 = qmAtomGroup[impropers[i].atom2]; 01666 int part3 = qmAtomGroup[impropers[i].atom3]; 01667 int part4 = qmAtomGroup[impropers[i].atom4]; 01668 if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0) { 01669 //CkPrintf("-----i %i IMPROPER ATOMS %i %i %i %i partitions %i %i %i %i\n",i,impropers[i].atom1, impropers[i].atom2, impropers[i].atom3, impropers[i].atom4, part1, part2, part3,part4); 01670 qmDroppedImpropers++; 01671 } 01672 else { 01673 nonQMImpropers[nonQMImproperCount++] = impropers[i]; 01674 } 01675 } 01676 numImpropers = nonQMImproperCount; 01677 delete [] impropers; 01678 impropers = new Improper[numImpropers]; 01679 for (int i = 0; i < numImpropers; i++) { 01680 impropers[i]=nonQMImpropers[i]; 01681 } 01682 delete [] nonQMImpropers; 01683 01684 // Crossterms 01685 Crossterm *nonQMCrossterms; 01686 nonQMCrossterms = new Crossterm[numCrossterms]; 01687 int nonQMCrosstermCount = 0; 01688 qmDroppedCrossterms = 0 ; 01689 for (int i = 0; i < numCrossterms; i++) { 01690 int part1 = qmAtomGroup[crossterms[i].atom1]; 01691 int part2 = qmAtomGroup[crossterms[i].atom2]; 01692 int part3 = qmAtomGroup[crossterms[i].atom3]; 01693 int part4 = qmAtomGroup[crossterms[i].atom4]; 01694 int part5 = qmAtomGroup[crossterms[i].atom5]; 01695 int part6 = qmAtomGroup[crossterms[i].atom6]; 01696 int part7 = qmAtomGroup[crossterms[i].atom7]; 01697 int part8 = qmAtomGroup[crossterms[i].atom8]; 01698 if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0 && 01699 part5 > 0 && part6 > 0 && part7 > 0 && part8 > 0) { 01700 //CkPrintf("-----i %i IMPROPER ATOMS %i %i %i %i partitions %i %i %i %i\n",i,impropers[i].atom1, impropers[i].atom2, impropers[i].atom3, impropers[i].atom4, part1, part2, part3,part4); 01701 qmDroppedCrossterms++; 01702 } 01703 else { 01704 nonQMCrossterms[nonQMCrosstermCount++] = crossterms[i]; 01705 } 01706 } 01707 numCrossterms = nonQMCrosstermCount; 01708 delete [] crossterms; 01709 crossterms = new Crossterm[numCrossterms] ; 01710 for (int i=0; i < numCrossterms; i++){ 01711 crossterms[i] = nonQMCrossterms[i] ; 01712 } 01713 delete [] nonQMCrossterms ; 01714 01715 if (!CkMyPe()) { 01716 iout << iINFO << "The QM region will remove " << qmDroppedBonds << " bonds, " << 01717 qmDroppedAngles << " angles, " << qmDroppedDihedrals << " dihedrals, " 01718 << qmDroppedImpropers << " impropers and " << qmDroppedCrossterms 01719 << " crossterms.\n" << endi ; 01720 } 01721 01722 #endif 01723 }
| void Molecule::freeBFactorData | ( | ) | [inline] |
Definition at line 1021 of file Molecule.h.
Referenced by NamdState::loadStructure().
| void Molecule::freeOccupancyData | ( | ) | [inline] |
Definition at line 1017 of file Molecule.h.
Referenced by NamdState::loadStructure().
| Bond* Molecule::get_acceptor | ( | int | dnum | ) | const [inline] |
Definition at line 1092 of file Molecule.h.
| Angle* Molecule::get_angle | ( | int | anum | ) | const [inline] |
Definition at line 1055 of file Molecule.h.
Referenced by atoms_1to4(), and dumpbench().
| int32* Molecule::get_angles_for_atom | ( | int | anum | ) | [inline] |
Definition at line 1134 of file Molecule.h.
Referenced by atoms_1to4(), and dumpbench().
| int Molecule::get_atom_from_index_in_residue | ( | const char * | segid, | |
| int | resid, | |||
| int | index | |||
| ) | const |
Definition at line 158 of file Molecule.C.
References ResidueLookupElem::lookup(), and NAMD_die().
Referenced by GlobalMasterFreeEnergy::getAtomID(), and GlobalMasterEasy::getAtomID().
00159 { 00160 00161 if (atomNames == NULL || resLookup == NULL) 00162 { 00163 NAMD_die("Tried to find atom from name on node other than node 0"); 00164 } 00165 int i = 0; 00166 int end = 0; 00167 if ( resLookup->lookup(segid,resid,&i,&end) ) return -1; 00168 if ( index >= 0 && index < ( end - i ) ) return ( index + i ); 00169 return -1; 00170 }
| int Molecule::get_atom_from_name | ( | const char * | segid, | |
| int | resid, | |||
| const char * | aname | |||
| ) | const |
Definition at line 121 of file Molecule.C.
References atomNamePool, ResidueLookupElem::lookup(), and NAMD_die().
Referenced by colvarproxy_namd::check_atom_id(), GlobalMasterFreeEnergy::getAtomID(), and GlobalMasterEasy::getAtomID().
00122 { 00123 00124 if (atomNames == NULL || resLookup == NULL) 00125 { 00126 NAMD_die("Tried to find atom from name on node other than node 0"); 00127 } 00128 00129 int i = 0; 00130 int end = 0; 00131 if ( resLookup->lookup(segid,resid,&i,&end) ) return -1; 00132 for ( ; i < end; ++i ) { 00133 #ifdef MEM_OPT_VERSION 00134 Index idx = atomNames[i].atomnameIdx; 00135 if(!strcasecmp(aname, atomNamePool[idx])) return i; 00136 #else 00137 if ( ! strcasecmp(aname,atomNames[i].atomname) ) return i; 00138 #endif 00139 } 00140 return -1; 00141 }
| const char* Molecule::get_atomtype | ( | int | anum | ) | const [inline] |
Definition at line 1103 of file Molecule.h.
References atomTypePool, and NAMD_die().
01104 { 01105 if (atomNames == NULL) 01106 { 01107 NAMD_die("Tried to find atom type on node other than node 0"); 01108 } 01109 01110 #ifdef MEM_OPT_VERSION 01111 return atomTypePool[atomNames[anum].atomtypeIdx]; 01112 #else 01113 return(atomNames[anum].atomtype); 01114 #endif 01115 }
| Bond* Molecule::get_bond | ( | int | bnum | ) | const [inline] |
Definition at line 1052 of file Molecule.h.
Referenced by atoms_1to4(), and dumpbench().
| int32* Molecule::get_bonds_for_atom | ( | int | anum | ) | [inline] |
Definition at line 1132 of file Molecule.h.
Referenced by atoms_1to4(), and dumpbench().
| int Molecule::get_cluster | ( | int | anum | ) | const [inline] |
Definition at line 1010 of file Molecule.h.
Referenced by wrap_coor_int().
| int Molecule::get_clusterSize | ( | int | anum | ) | const [inline] |
Definition at line 1011 of file Molecule.h.
Referenced by wrap_coor_int().
Definition at line 1249 of file Molecule.h.
Referenced by ComputeRestraints::doForce().
| void Molecule::get_constorque_params | ( | BigReal & | v, | |
| Vector & | a, | |||
| Vector & | p, | |||
| int | atomnum | |||
| ) | const [inline] |
Definition at line 1323 of file Molecule.h.
References consTorqueIndexes, and consTorqueParams.
Referenced by ComputeConsTorque::doForce().
01325 { 01326 v = consTorqueParams[consTorqueIndexes[atomnum]].v; 01327 a = consTorqueParams[consTorqueIndexes[atomnum]].a; 01328 p = consTorqueParams[consTorqueIndexes[atomnum]].p; 01329 }
| Crossterm* Molecule::get_crossterm | ( | int | inum | ) | const [inline] |
Definition at line 1064 of file Molecule.h.
| int32* Molecule::get_crossterms_for_atom | ( | int | anum | ) | [inline] |
Definition at line 1140 of file Molecule.h.
| const Real* Molecule::get_cSMDcoffs | ( | ) | [inline] |
Definition at line 866 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| int const* const* Molecule::get_cSMDindex | ( | ) | [inline] |
Definition at line 862 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| int const* Molecule::get_cSMDindxLen | ( | ) | [inline] |
Definition at line 861 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| const Real* Molecule::get_cSMDKs | ( | ) | [inline] |
Definition at line 864 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| int Molecule::get_cSMDnumInst | ( | ) | [inline] |
Definition at line 860 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| int const* const* Molecule::get_cSMDpairs | ( | ) | [inline] |
Definition at line 863 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| const Real* Molecule::get_cSMDVels | ( | ) | [inline] |
Definition at line 865 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| Dihedral* Molecule::get_dihedral | ( | int | dnum | ) | const [inline] |
Definition at line 1061 of file Molecule.h.
Referenced by atoms_1to4(), and dumpbench().
| int32* Molecule::get_dihedrals_for_atom | ( | int | anum | ) | [inline] |
Definition at line 1136 of file Molecule.h.
Referenced by atoms_1to4(), and dumpbench().
| Bond* Molecule::get_donor | ( | int | dnum | ) | const [inline] |
Definition at line 1089 of file Molecule.h.
| const ExclusionCheck* Molecule::get_excl_check_for_atom | ( | int | anum | ) | const [inline] |
Definition at line 1160 of file Molecule.h.
Referenced by dumpbench().
| Exclusion* Molecule::get_exclusion | ( | int | ex | ) | const [inline] |
Definition at line 1099 of file Molecule.h.
References exclusions.
| int32* Molecule::get_exclusions_for_atom | ( | int | anum | ) | [inline] |
Definition at line 1142 of file Molecule.h.
References exclusionsByAtom.
| int Molecule::get_fep_bonded_type | ( | const int * | atomID, | |
| unsigned int | order | |||
| ) | const [inline] |
Definition at line 1372 of file Molecule.h.
References NAMD_die(), and simParams.
Referenced by ImproperElem::computeForce(), DihedralElem::computeForce(), BondElem::computeForce(), AnisoElem::computeForce(), and AngleElem::computeForce().
01373 { 01374 int typeSum = 0; 01375 for ( int i=0; i < order; ++i ) { 01376 typeSum += (fepAtomFlags[atomID[i]] == 2 ? -1 : fepAtomFlags[atomID[i]]); 01377 } 01378 // Increase the cutoff if scaling purely alchemical bonds. 01379 // This really only applies when order = 2. 01380 if ( simParams->alchBondDecouple ) order++; 01381 01382 // The majority of interactions are type 0, so catch those first. 01383 if ( typeSum == 0 || abs(typeSum) == order ) return 0; 01384 else if ( 0 < typeSum && typeSum < order ) return 1; 01385 else if ( -order < typeSum && typeSum < 0 ) return 2; 01386 01387 // Alchemify should always keep this from bombing, but just in case... 01388 NAMD_die("Unexpected alchemical bonded interaction!"); 01389 return 0; 01390 }
| unsigned char Molecule::get_fep_type | ( | int | anum | ) | const [inline] |
Definition at line 1332 of file Molecule.h.
Referenced by TholeElem::computeForce(), CrosstermElem::computeForce(), WorkDistrib::createAtomLists(), and ComputeHomeTuples< TholeElem, Thole, TholeValue >::loadTuples().
| const int32* Molecule::get_full_exclusions_for_atom | ( | int | anum | ) | const [inline] |
Definition at line 1144 of file Molecule.h.
Referenced by ComputeNonbondedCUDA::build_exclusions().
| Real Molecule::get_go_cutoff | ( | int | chain1, | |
| int | chain2 | |||
| ) | [inline] |
Definition at line 1588 of file Molecule.h.
References go_val::cutoff, go_array, and MAX_GO_CHAINS.
Referenced by build_go_arrays(), build_go_sigmas(), build_go_sigmas2(), get_go_force(), get_go_force2(), and get_go_force_new().
01588 { return go_array[MAX_GO_CHAINS*chain1 + chain2].cutoff; };
| Real Molecule::get_go_epsilon | ( | int | chain1, | |
| int | chain2 | |||
| ) | [inline] |
Definition at line 1597 of file Molecule.h.
References go_val::epsilon, go_array, and MAX_GO_CHAINS.
Referenced by get_go_force(), get_go_force2(), and get_go_force_new().
01597 { return go_array[MAX_GO_CHAINS*chain1 + chain2].epsilon; };
| Real Molecule::get_go_epsilonRep | ( | int | chain1, | |
| int | chain2 | |||
| ) | [inline] |
Definition at line 1591 of file Molecule.h.
References go_val::epsilonRep, go_array, and MAX_GO_CHAINS.
Referenced by get_go_force(), get_go_force2(), and get_go_force_new().
01591 { return go_array[MAX_GO_CHAINS*chain1 + chain2].epsilonRep; };
| int Molecule::get_go_exp_a | ( | int | chain1, | |
| int | chain2 | |||
| ) | [inline] |
Definition at line 1600 of file Molecule.h.
References go_val::exp_a, go_array, and MAX_GO_CHAINS.
Referenced by build_go_sigmas(), build_go_sigmas2(), get_go_force(), get_go_force2(), and get_go_force_new().
01600 { return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_a; };
| int Molecule::get_go_exp_b | ( | int | chain1, | |
| int | chain2 | |||
| ) | [inline] |
Definition at line 1603 of file Molecule.h.
References go_val::exp_b, go_array, and MAX_GO_CHAINS.
Referenced by build_go_sigmas(), build_go_sigmas2(), get_go_force(), get_go_force2(), and get_go_force_new().
01603 { return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_b; };
| int Molecule::get_go_exp_rep | ( | int | chain1, | |
| int | chain2 | |||
| ) | [inline] |
Definition at line 1606 of file Molecule.h.
References go_val::exp_rep, go_array, and MAX_GO_CHAINS.
Referenced by get_go_force(), get_go_force2(), and get_go_force_new().
01606 { return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_rep; };
| BigReal Molecule::get_go_force | ( | BigReal | r, | |
| int | atom1, | |||
| int | atom2, | |||
| BigReal * | goNative, | |||
| BigReal * | goNonnative | |||
| ) | const |
Definition at line 1260 of file GoMolecule.C.
References atomChainTypes, get_go_cutoff(), get_go_epsilon(), get_go_epsilonRep(), get_go_exp_a(), get_go_exp_b(), get_go_exp_rep(), get_go_sigmaRep(), goSigmaIndices, goSigmas, goWithinCutoff, and numGoAtoms.
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 }
| BigReal Molecule::get_go_force2 | ( | BigReal | x, | |
| BigReal | y, | |||
| BigReal | z, | |||
| int | atom1, | |||
| int | atom2, | |||
| BigReal * | goNative, | |||
| BigReal * | goNonnative | |||
| ) | const |
Definition at line 1456 of file GoMolecule.C.
References atomChainTypes, get_go_cutoff(), get_go_epsilon(), get_go_epsilonRep(), get_go_exp_a(), get_go_exp_b(), get_go_exp_rep(), get_go_sigmaRep(), go_restricted(), goIndxLJB, goResidIndices, goSigmaPairA, goSigmaPairB, pointerToGoBeg, and pointerToGoEnd.
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 }
| BigReal Molecule::get_go_force_new | ( | BigReal | r, | |
| int | atom1, | |||
| int | atom2, | |||
| BigReal * | goNative, | |||
| BigReal * | goNonnative | |||
| ) | const |
Definition at line 1334 of file GoMolecule.C.
References atomChainTypes, DebugM, get_go_cutoff(), get_go_epsilon(), get_go_epsilonRep(), get_go_exp_a(), get_go_exp_b(), get_go_exp_rep(), get_go_sigmaRep(), go_restricted(), goCoordinates, goResids, and goSigmaIndices.
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/(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 //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 }
| Real Molecule::get_go_sigmaRep | ( | int | chain1, | |
| int | chain2 | |||
| ) | [inline] |
Definition at line 1594 of file Molecule.h.
References go_array, MAX_GO_CHAINS, and go_val::sigmaRep.
Referenced by get_go_force(), get_go_force2(), and get_go_force_new().
01594 { return go_array[MAX_GO_CHAINS*chain1 + chain2].sigmaRep; };
| GridforceGrid* Molecule::get_gridfrc_grid | ( | int | gridnum | ) | const [inline] |
Definition at line 1263 of file Molecule.h.
References numGridforceGrids.
Referenced by ComputeGridForce::doForce(), Node::reloadGridforceGrid(), and Node::updateGridScale().
01264 { 01265 GridforceGrid *result = NULL; 01266 if (gridnum >= 0 && gridnum < numGridforceGrids) { 01267 result = gridfrcGrid[gridnum]; 01268 } 01269 return result; 01270 }
Definition at line 1257 of file Molecule.h.
Referenced by ComputeGridForce::do_calc().
| BigReal Molecule::get_gro_force2 | ( | BigReal | x, | |
| BigReal | y, | |||
| BigReal | z, | |||
| int | atom1, | |||
| int | atom2, | |||
| BigReal * | pairLJEnergy, | |||
| BigReal * | pairGaussEnergy | |||
| ) | const |
Definition at line 1173 of file GoMolecule.C.
References A, gA, giSigma1, giSigma2, gMu1, gMu2, gRepulsive, indxGaussB, indxLJB, pairC12, pairC6, pointerToGaussBeg, pointerToGaussEnd, pointerToLJBeg, and pointerToLJEnd.
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 }
| int Molecule::get_groupSize | ( | int | ) |
| Improper* Molecule::get_improper | ( | int | inum | ) | const [inline] |
Definition at line 1058 of file Molecule.h.
Referenced by dumpbench().
| int32* Molecule::get_impropers_for_atom | ( | int | anum | ) | [inline] |
Definition at line 1138 of file Molecule.h.
Referenced by dumpbench().
| Lphost* Molecule::get_lphost | ( | int | atomid | ) | const [inline] |
Definition at line 1068 of file Molecule.h.
| const int32* Molecule::get_mod_exclusions_for_atom | ( | int | anum | ) | const [inline] |
Definition at line 1146 of file Molecule.h.
| int Molecule::get_mother_atom | ( | int | ) |
| void Molecule::get_movdrag_params | ( | Vector & | v, | |
| int | atomnum | |||
| ) | const [inline] |
Definition at line 1308 of file Molecule.h.
Referenced by Sequencer::addMovDragToPosition().
| Bool Molecule::get_noPC | ( | ) | [inline] |
Definition at line 848 of file Molecule.h.
Referenced by ComputeQM::initialize().
| int Molecule::get_numQMAtoms | ( | ) | [inline] |
Definition at line 817 of file Molecule.h.
Referenced by ComputeQMMgr::calcMOPAC(), ComputeQMMgr::calcORCA(), ComputeQMMgr::calcUSR(), ComputePme::doQMWork(), ComputeQM::initialize(), ComputeQMMgr::procQMRes(), HomePatch::qmSwapAtoms(), ComputeQMMgr::recvPartQM(), and ComputeQM::saveResults().
| Real* Molecule::get_qmAtmChrg | ( | ) | [inline] |
Definition at line 814 of file Molecule.h.
Referenced by ComputeQMMgr::calcMOPAC(), ComputeQMMgr::calcORCA(), ComputeQMMgr::calcUSR(), ComputePme::doQMWork(), ComputeQM::initialize(), HomePatch::qmSwapAtoms(), and ComputeQM::saveResults().
| const int* Molecule::get_qmAtmIndx | ( | ) | [inline] |
Definition at line 815 of file Molecule.h.
Referenced by ComputeQMMgr::calcMOPAC(), ComputeQMMgr::calcORCA(), ComputeQMMgr::calcUSR(), ComputePme::doQMWork(), ComputeQM::initialize(), HomePatch::qmSwapAtoms(), and ComputeQM::saveResults().
| Real Molecule::get_qmAtomGroup | ( | int | indx | ) | const [inline] |
Definition at line 812 of file Molecule.h.
| const Real* Molecule::get_qmAtomGroup | ( | ) | const [inline] |
Definition at line 811 of file Molecule.h.
Referenced by ComputeQMMgr::calcMOPAC(), ComputeQMMgr::calcORCA(), ComputeQMMgr::calcUSR(), ComputePme::doQMWork(), ComputeQM::initialize(), HomePatch::qmSwapAtoms(), and ComputeQM::saveResults().
| Bool Molecule::get_qmcSMD | ( | ) | [inline] |
Definition at line 859 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| int* Molecule::get_qmCustomPCIdxs | ( | ) | [inline] |
Definition at line 857 of file Molecule.h.
Referenced by ComputeQM::initialize().
| int* Molecule::get_qmCustPCSizes | ( | ) | [inline] |
Definition at line 856 of file Molecule.h.
Referenced by ComputeQM::initialize().
| const BigReal* Molecule::get_qmDummyBondVal | ( | ) | [inline] |
Definition at line 826 of file Molecule.h.
Referenced by ComputeQMMgr::recvPntChrg().
| const char* const* Molecule::get_qmDummyElement | ( | ) | [inline] |
Definition at line 831 of file Molecule.h.
Referenced by ComputeQMMgr::recvPntChrg().
| const char* const* Molecule::get_qmElements | ( | ) | [inline] |
Definition at line 818 of file Molecule.h.
Referenced by ComputeQMMgr::recvPntChrg().
| const int* const* Molecule::get_qmGrpBonds | ( | ) | [inline] |
Definition at line 829 of file Molecule.h.
Referenced by ComputeQMMgr::recvPntChrg().
| Real* Molecule::get_qmGrpChrg | ( | ) | [inline] |
Definition at line 822 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| Real* Molecule::get_qmGrpID | ( | ) | [inline] |
Definition at line 821 of file Molecule.h.
Referenced by ComputeQM::initialize(), and ComputeQMMgr::recvPartQM().
| Real* Molecule::get_qmGrpMult | ( | ) | [inline] |
Definition at line 823 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| const int* Molecule::get_qmGrpNumBonds | ( | ) | [inline] |
Definition at line 827 of file Molecule.h.
Referenced by ComputeQMMgr::recvPntChrg().
| const int* Molecule::get_qmGrpSizes | ( | ) | [inline] |
Definition at line 820 of file Molecule.h.
Referenced by ComputeQMMgr::recvPntChrg().
| int Molecule::get_qmLSSFreq | ( | ) | [inline] |
Definition at line 839 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| int* Molecule::get_qmLSSIdxs | ( | ) | [inline] |
Definition at line 837 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| Mass* Molecule::get_qmLSSMass | ( | ) | [inline] |
Definition at line 838 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| int* Molecule::get_qmLSSRefIDs | ( | ) | [inline] |
Definition at line 841 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| Mass* Molecule::get_qmLSSRefMass | ( | ) | [inline] |
Definition at line 843 of file Molecule.h.
| int* Molecule::get_qmLSSRefSize | ( | ) | [inline] |
Definition at line 842 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| int Molecule::get_qmLSSResSize | ( | ) | [inline] |
Definition at line 840 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| int* Molecule::get_qmLSSSize | ( | ) | [inline] |
Definition at line 836 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| int* Molecule::get_qmMeMMindx | ( | ) | [inline] |
Definition at line 850 of file Molecule.h.
Referenced by ComputeQM::initialize().
| int Molecule::get_qmMeNumBonds | ( | ) | [inline] |
Definition at line 849 of file Molecule.h.
Referenced by ComputeQM::initialize(), and ComputeQMMgr::recvPartQM().
| Real* Molecule::get_qmMeQMGrp | ( | ) | [inline] |
Definition at line 851 of file Molecule.h.
Referenced by ComputeQM::initialize().
| const int* const* Molecule::get_qmMMBond | ( | ) | [inline] |
Definition at line 828 of file Molecule.h.
Referenced by ComputeQMMgr::recvPntChrg().
| const int* const* Molecule::get_qmMMBondedIndx | ( | ) | [inline] |
Definition at line 830 of file Molecule.h.
Referenced by ComputeQMMgr::recvPntChrg().
| const int* const* Molecule::get_qmMMChargeTarget | ( | ) | [inline] |
Definition at line 833 of file Molecule.h.
Referenced by ComputeQMMgr::recvPntChrg().
| const int* Molecule::get_qmMMNumTargs | ( | ) | [inline] |
Definition at line 834 of file Molecule.h.
Referenced by ComputeQMMgr::recvPntChrg().
| std::map<int,int>& Molecule::get_qmMMSolv | ( | ) | [inline] |
Definition at line 844 of file Molecule.h.
| int Molecule::get_qmNumBonds | ( | ) | [inline] |
Definition at line 825 of file Molecule.h.
| int Molecule::get_qmNumGrps | ( | ) | const [inline] |
Definition at line 819 of file Molecule.h.
Referenced by ComputeQM::initialize(), ComputeQMMgr::recvPartQM(), and ComputeQM::saveResults().
| int Molecule::get_qmPCFreq | ( | ) | [inline] |
Definition at line 853 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| Bool Molecule::get_qmReplaceAll | ( | ) | [inline] |
Definition at line 846 of file Molecule.h.
Referenced by ComputeQMMgr::recvPartQM().
| int Molecule::get_qmTotCustPCs | ( | ) | [inline] |
Definition at line 855 of file Molecule.h.
Referenced by ComputeQM::initialize().
| int Molecule::get_residue_size | ( | const char * | segid, | |
| int | resid | |||
| ) | const |
Definition at line 144 of file Molecule.C.
References ResidueLookupElem::lookup(), and NAMD_die().
Referenced by GlobalMasterFreeEnergy::getNumAtoms(), GlobalMasterEasy::getNumAtoms(), and prepare_qm().
| void Molecule::get_rotdrag_params | ( | BigReal & | v, | |
| Vector & | a, | |||
| Vector & | p, | |||
| int | atomnum | |||
| ) | const [inline] |
Definition at line 1314 of file Molecule.h.
Referenced by Sequencer::addRotDragToPosition().
| unsigned char Molecule::get_ss_type | ( | int | anum | ) | const [inline] |
Definition at line 1338 of file Molecule.h.
Referenced by ComputeHomeTuples< TholeElem, Thole, TholeValue >::loadTuples(), and Sequencer::rescaleSoluteCharges().
| void Molecule::get_stir_refPos | ( | Vector & | refPos, | |
| int | atomnum | |||
| ) | const [inline] |
Definition at line 1289 of file Molecule.h.
| Real Molecule::get_stir_startTheta | ( | int | atomnum | ) | const [inline] |
Definition at line 1301 of file Molecule.h.
Referenced by ComputeStir::doForce().
| Bond* Molecule::getAllAcceptors | ( | ) | const [inline] |
Definition at line 1095 of file Molecule.h.
| Angle* Molecule::getAllAngles | ( | ) | const [inline] |
| Bond* Molecule::getAllBonds | ( | ) | const [inline] |
| Crossterm* Molecule::getAllCrossterms | ( | ) | const [inline] |
Definition at line 1081 of file Molecule.h.
Referenced by buildCrosstermData().
| Dihedral* Molecule::getAllDihedrals | ( | ) | const [inline] |
Definition at line 1080 of file Molecule.h.
Referenced by buildDihedralData().
| Bond* Molecule::getAllDonors | ( | ) | const [inline] |
Definition at line 1094 of file Molecule.h.
| Improper* Molecule::getAllImpropers | ( | ) | const [inline] |
Definition at line 1079 of file Molecule.h.
Referenced by buildImproperData().
| Lphost* Molecule::getAllLphosts | ( | ) | const [inline] |
Definition at line 1085 of file Molecule.h.
| BigReal Molecule::GetAtomAlpha | ( | int | i | ) | const [inline] |
Definition at line 480 of file Molecule.h.
References drude_constants::alpha.
00480 { 00481 return drudeConsts[i].alpha; 00482 }
| AtomNameInfo* Molecule::getAtomNames | ( | ) | const [inline] |
Definition at line 489 of file Molecule.h.
Referenced by buildAtomData().
| Atom* Molecule::getAtoms | ( | ) | const [inline] |
Definition at line 488 of file Molecule.h.
References atoms.
Referenced by buildAtomData(), WorkDistrib::createAtomLists(), and outputCompressedFile().
| AtomSegResInfo* Molecule::getAtomSegResInfo | ( | ) | const [inline] |
Definition at line 492 of file Molecule.h.
Referenced by buildAtomData().
| const float* Molecule::getBFactorData | ( | ) | [inline] |
Definition at line 1019 of file Molecule.h.
Referenced by NamdState::loadStructure(), and outputCompressedFile().
| BigReal Molecule::getEnergyTailCorr | ( | const | BigReal, | |
| const | int | |||
| ) |
Referenced by Controller::printEnergies().
| int const* Molecule::getLcpoParamType | ( | ) | [inline] |
Definition at line 476 of file Molecule.h.
Referenced by HomePatch::setLcpoType().
| const float* Molecule::getOccupancyData | ( | ) | [inline] |
Definition at line 1015 of file Molecule.h.
Referenced by NamdState::loadStructure(), and outputCompressedFile().
| BigReal Molecule::getVirialTailCorr | ( | const | BigReal | ) |
Referenced by Controller::calcPressure().
| Bool Molecule::go_restricted | ( | int | chain1, | |
| int | chain2, | |||
| int | rDiff | |||
| ) |
Definition at line 525 of file GoMolecule.C.
References FALSE, go_array, MAX_GO_CHAINS, MAX_RESTRICTIONS, and TRUE.
Referenced by build_go_arrays(), build_go_sigmas(), build_go_sigmas2(), get_go_force2(), and get_go_force_new().
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 }
| void Molecule::goInit | ( | ) |
Definition at line 54 of file GoMolecule.C.
References atomChainTypes, energyNative, energyNonnative, goCoordinates, goPDB, goResids, goSigmaIndices, goSigmas, goWithinCutoff, and numGoAtoms.
00054 { 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 }
| void Molecule::initialize | ( | ) |
Referenced by Molecule().
| Bool Molecule::is_atom_constorqued | ( | int | atomnum | ) | const [inline] |
Definition at line 1233 of file Molecule.h.
References consTorqueIndexes, FALSE, and numConsTorque.
01234 { 01235 if (numConsTorque) 01236 { 01237 // Check the index to see if it is constrained 01238 return(consTorqueIndexes[atomnum] != -1); 01239 } 01240 else 01241 { 01242 // No constraints at all, so just return FALSE 01243 return(FALSE); 01244 } 01245 }
| Bool Molecule::is_atom_constrained | ( | int | atomnum | ) | const [inline] |
Definition at line 1184 of file Molecule.h.
References FALSE, and numConstraints.
Referenced by ComputeRestraints::doForce().
01185 { 01186 if (numConstraints) 01187 { 01188 // Check the index to see if it is constrained 01189 return(consIndexes[atomnum] != -1); 01190 } 01191 else 01192 { 01193 // No constraints at all, so just return FALSE 01194 return(FALSE); 01195 } 01196 }
| Bool Molecule::is_atom_exPressure | ( | int | atomnum | ) | const [inline] |
Definition at line 1434 of file Molecule.h.
References numExPressureAtoms.
Referenced by Sequencer::langevinPiston().
01435 { 01436 return (numExPressureAtoms && exPressureAtomFlags[atomnum]); 01437 }
| Bool Molecule::is_atom_fixed | ( | int | atomnum | ) | const [inline] |
Definition at line 1394 of file Molecule.h.
References numFixedAtoms.
Referenced by WorkDistrib::createAtomLists().
01395 { 01396 return (numFixedAtoms && fixedAtomFlags[atomnum]); 01397 }
| Bool Molecule::is_atom_gridforced | ( | int | atomnum, | |
| int | gridnum | |||
| ) | const [inline] |
Definition at line 1168 of file Molecule.h.
References FALSE, and numGridforceGrids.
Referenced by ComputeGridForce::do_calc().
01169 { 01170 if (numGridforceGrids) 01171 { 01172 return(gridfrcIndexes[gridnum][atomnum] != -1); 01173 } 01174 else 01175 { 01176 return(FALSE); 01177 } 01178 }
| Bool Molecule::is_atom_movdragged | ( | int | atomnum | ) | const [inline] |
Definition at line 1201 of file Molecule.h.
References FALSE, and numMovDrag.
Referenced by Sequencer::addMovDragToPosition().
01202 { 01203 if (numMovDrag) 01204 { 01205 // Check the index to see if it is constrained 01206 return(movDragIndexes[atomnum] != -1); 01207 } 01208 else 01209 { 01210 // No constraints at all, so just return FALSE 01211 return(FALSE); 01212 } 01213 }
| Bool Molecule::is_atom_rotdragged | ( | int | atomnum | ) | const [inline] |
Definition at line 1217 of file Molecule.h.
References FALSE, and numRotDrag.
Referenced by Sequencer::addRotDragToPosition().
01218 { 01219 if (numRotDrag) 01220 { 01221 // Check the index to see if it is constrained 01222 return(rotDragIndexes[atomnum] != -1); 01223 } 01224 else 01225 { 01226 // No constraints at all, so just return FALSE 01227 return(FALSE); 01228 } 01229 }
| Bool Molecule::is_atom_stirred | ( | int | atomnum | ) | const [inline] |
Definition at line 1415 of file Molecule.h.
References FALSE, and numStirredAtoms.
Referenced by ComputeStir::doForce().
01416 { 01417 if (numStirredAtoms) 01418 { 01419 // Check the index to see if it is constrained 01420 return(stirIndexes[atomnum] != -1); 01421 } 01422 else 01423 { 01424 // No constraints at all, so just return FALSE 01425 return(FALSE); 01426 } 01427 }
| Bool Molecule::is_drude | ( | int | ) |
| Bool Molecule::is_group_fixed | ( | int | atomnum | ) | const [inline] |
Definition at line 1430 of file Molecule.h.
References numFixedAtoms.
01431 { 01432 return (numFixedAtoms && (fixedAtomFlags[atomnum] == -1)); 01433 }
| Bool Molecule::is_hydrogen | ( | int | ) |
| Bool Molecule::is_hydrogenGroupParent | ( | int | ) |
| Bool Molecule::is_lp | ( | int | ) |
| Bool Molecule::is_oxygen | ( | int | ) |
| Bool Molecule::is_water | ( | int | ) |
Referenced by outputCompressedFile(), and wrap_coor_int().
| Real Molecule::langevin_param | ( | int | atomnum | ) | const [inline] |
Definition at line 1283 of file Molecule.h.
Referenced by WorkDistrib::createAtomLists().
| int64_t Molecule::num_deg_freedom | ( | int | isInitialReport = 0 |
) | const [inline] |
Definition at line 521 of file Molecule.h.
References num_fixed_atoms(), numAtoms, numConstraints, numFepInitial, numFixedRigidBonds, numLonepairs, numRigidBonds, simParams, and WAT_TIP4.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), Controller::Controller(), NamdState::loadStructure(), and Controller::receivePressure().
00521 { 00522 // local variables prefixed by s_ 00523 int64_t s_NumDegFreedom = 3 * (int64_t) numAtoms; 00524 int64_t s_NumFixedAtoms = num_fixed_atoms(); 00525 if (s_NumFixedAtoms) s_NumDegFreedom -= 3 * s_NumFixedAtoms; 00526 if (numLonepairs) s_NumDegFreedom -= 3 * numLonepairs; 00527 if ( ! (s_NumFixedAtoms || numConstraints 00528 || simParams->comMove || simParams->langevinOn) ) { 00529 s_NumDegFreedom -= 3; 00530 } 00531 if ( ! isInitialReport && simParams->pairInteractionOn) { 00532 // 00533 // DJH: a kludge? We want to initially report system degrees of freedom 00534 // 00535 // this doesn't attempt to deal with fixed atoms or constraints 00536 s_NumDegFreedom = 3 * numFepInitial; 00537 } 00538 int s_NumFixedRigidBonds = 00539 (simParams->fixedAtomsOn ? numFixedRigidBonds : 0); 00540 if (simParams->watmodel == WAT_TIP4) { 00541 // numLonepairs is subtracted here because all lonepairs have a rigid bond 00542 // to oxygen, but all of the LP degrees of freedom are dealt with above 00543 s_NumDegFreedom -= (numRigidBonds - s_NumFixedRigidBonds - numLonepairs); 00544 } 00545 else { 00546 // Non-polarized systems don't have LPs. 00547 // For Drude model, bonds that attach LPs are not counted as rigid; 00548 // LPs have already been subtracted from degrees of freedom. 00549 s_NumDegFreedom -= (numRigidBonds - s_NumFixedRigidBonds); 00550 } 00551 return s_NumDegFreedom; 00552 }
| int Molecule::num_fixed_atoms | ( | ) | const [inline] |
Definition at line 495 of file Molecule.h.
References numFixedAtoms, and simParams.
Referenced by num_deg_freedom(), num_fixed_groups(), num_group_deg_freedom(), and Controller::receivePressure().
00495 { 00496 // local variables prefixed by s_ 00497 int s_NumFixedAtoms = (simParams->fixedAtomsOn ? numFixedAtoms : 0); 00498 return s_NumFixedAtoms; // value is "turned on" SimParameters 00499 }
| int Molecule::num_fixed_groups | ( | ) | const [inline] |
Definition at line 501 of file Molecule.h.
References num_fixed_atoms(), and numFixedGroups.
Referenced by num_group_deg_freedom(), and Controller::receivePressure().
00501 { 00502 // local variables prefixed by s_ 00503 int s_NumFixedAtoms = num_fixed_atoms(); 00504 int s_NumFixedGroups = (s_NumFixedAtoms ? numFixedGroups : 0); 00505 return s_NumFixedGroups; // value is "turned on" SimParameters 00506 }
| int64_t Molecule::num_group_deg_freedom | ( | ) | const [inline] |
Definition at line 508 of file Molecule.h.
References num_fixed_atoms(), num_fixed_groups(), numConstraints, numHydrogenGroups, and simParams.
Referenced by Controller::receivePressure().
00508 { 00509 // local variables prefixed by s_ 00510 int64_t s_NumGroupDegFreedom = 3 * (int64_t) numHydrogenGroups; 00511 int64_t s_NumFixedAtoms = num_fixed_atoms(); 00512 int s_NumFixedGroups = num_fixed_groups(); 00513 if (s_NumFixedGroups) s_NumGroupDegFreedom -= 3 * s_NumFixedGroups; 00514 if ( ! (s_NumFixedAtoms || numConstraints 00515 || simParams->comMove || simParams->langevinOn) ) { 00516 s_NumGroupDegFreedom -= 3; 00517 } 00518 return s_NumGroupDegFreedom; 00519 }
| void Molecule::prepare_qm | ( | const char * | pdbFileName, | |
| Parameters * | params, | |||
| ConfigList * | cfgList | |||
| ) |
Prepares Live Solvent Selection
Data gathering from PDB to find QM atom and bond info
Multiplicity of each QM group
Charge of each QM group
Populate arrays that are used throughout the the calculations.
Overides Link Atom element with user selection.
Bond Schemes. Prepares for treatment of QM-MM bonds in ComputeQM.C
Live Solvent Selection
Custom Point Charge selection
Topology preparation
Definition at line 109 of file MoleculeQM.C.
References ResizeArray< Elem >::add(), Parameters::assign_vdw_index(), qmSolvData::atmIDs, PDB::atom(), bond::atom1, bond::atom2, ResizeArray< Elem >::begin(), bond::bond_type, atom_constants::charge, StringList::data, DebugM, PDBAtom::element(), ResizeArray< Elem >::end(), endi(), SimParameters::extraBondsOn, SortedArray< Elem >::find(), ConfigList::find(), SimParameters::fixedAtomsOn, get_residue_size(), ObjectArena< Type >::getNewArray(), iERROR(), iINFO(), SortedArray< Elem >::insert(), iout, iWARN(), ResidueLookupElem::lookup(), atom_constants::mass, NAMD_blank_string(), NAMD_die(), NAMD_read_line(), StringList::next, PDB::num_atoms(), numAtoms, numBonds, numRealBonds, PDBAtom::occupancy(), SimParameters::PMEOn, SimParameters::qmBondDist, SimParameters::qmBondOn, SimParameters::qmBondScheme, SimParameters::qmChrgFromPSF, SimParameters::qmColumn, SimParameters::qmCSMD, SimParameters::qmCustomPCSel, SimParameters::qmElecEmbed, SimParameters::qmFormat, QMFormatMOPAC, QMFormatORCA, SimParameters::qmLSSFreq, SimParameters::qmLSSMode, QMLSSMODECOM, QMLSSMODEDIST, SimParameters::qmLSSOn, SimParameters::qmLSSResname, SimParameters::qmMOPACAddConfigChrg, SimParameters::qmNoPC, SimParameters::qmPCSelFreq, QMSCHEMECS, QMSCHEMERCD, QMSCHEMEZ1, QMSCHEMEZ2, QMSCHEMEZ3, SimParameters::qmVDW, PDBAtom::residuename(), PDBAtom::residueseq(), PDBAtom::segmentname(), ResizeArray< Elem >::size(), split(), SimParameters::stepsPerCycle, and PDBAtom::temperaturefactor().
Referenced by NamdState::loadStructure().
00110 { 00111 #ifdef MEM_OPT_VERSION 00112 NAMD_die("QMMM interface is not supported in memory optimized builds."); 00113 #else 00114 00115 PDB *pdbP; // Pointer to PDB object to use 00116 00117 qmNumQMAtoms = 0; 00118 00119 qmNumGrps = 0 ; 00120 00121 iout << iINFO << "Using the following PDB file for QM parameters: " << 00122 pdbFileName << "\n" << endi; 00123 if (pdbFileName) 00124 pdbP = new PDB(pdbFileName); 00125 else 00126 iout << "pdbFile name not defined!\n\n" << endi ; 00127 00128 if (pdbP->num_atoms() != numAtoms) 00129 { 00130 NAMD_die("Number of atoms in QM parameter PDB file doesn't match coordinate PDB"); 00131 } 00132 00133 qmElementArray = new char*[numAtoms] ; 00134 00135 qmAtomGroup = new Real[numAtoms] ; 00136 00137 BigReal bondVal; 00138 Real atmGrp; 00139 00140 std::set<Real> atmGrpSet; 00141 00142 std::vector<Real> grpChrgVec; 00143 00144 // Value in the bond column fro QM-MM bonds. 00145 std::vector<BigReal> qmBondValVec; 00146 // Identifiers of different QM regions 00147 std::vector<Real> qmGroupIDsVec; 00148 // This maps the qm group ID with the serail index of this group in 00149 // various arrays. 00150 std::map<Real,int> qmGrpIDMap ; 00151 00152 std::vector<int> qmAtmIndxVec, qmGrpSizeVec; 00153 00154 // Used to parse user selection in different places. 00155 std::vector<std::string> strVec; 00156 00157 qmNoPC = simParams->qmNoPC; 00158 00159 // We set to zero by default, So there is no need for extra processing. 00160 qmPCFreq = 0; 00161 if (simParams->qmPCSelFreq > 1) 00162 qmPCFreq = simParams->qmPCSelFreq; 00163 00164 00167 00168 00169 qmLSSTotalNumAtms = 0; 00170 SortedArray< qmSolvData> lssGrpRes; 00171 std::map<Real,std::vector<int> > lssGrpRefIDs; 00172 refSelStrMap lssRefUsrSel; 00173 int totalNumRefAtms = 0; 00174 int lssClassicResIndx = 0 ; 00175 int lssCurrClassResID = -1 ; 00176 char lssCurrClassResSegID[5]; 00177 if (simParams->qmLSSOn) { 00178 DebugM(4, "LSS is ON! Processing QM solvent.\n") ; 00179 00180 if (resLookup == NULL) 00181 NAMD_die("We need residue data to conduct LSS."); 00182 00183 if (simParams->qmLSSMode == QMLSSMODECOM) { 00184 00185 StringList *current = cfgList->find("QMLSSRef"); 00186 for ( ; current; current = current->next ) { 00187 00188 strVec = split( std::string(current->data) , " "); 00189 00190 if (strVec.size() != 3 ) { 00191 iout << iERROR << "Format error in QM LSS size: " 00192 << current->data 00193 << "\n" << endi; 00194 NAMD_die("Error processing QM information."); 00195 } 00196 00197 std::stringstream storConv ; 00198 00199 storConv << strVec[0] ; 00200 Real grpID ; 00201 storConv >> grpID; 00202 if (storConv.fail()) { 00203 iout << iERROR << "Error parsing QMLSSRef selection: " 00204 << current->data 00205 << "\n" << endi; 00206 NAMD_die("Error processing QM information."); 00207 } 00208 00209 std::string segName = strVec[1].substr(0,4); 00210 00211 storConv.clear() ; 00212 storConv << strVec[2]; 00213 int resID ; 00214 storConv >> resID; 00215 if (storConv.fail()) { 00216 iout << iERROR << "Error parsing QMLSSRef selection: " 00217 << current->data 00218 << "\n" << endi; 00219 NAMD_die("Error processing QM information."); 00220 } 00221 00222 auto it = lssRefUsrSel.find(grpID) ; 00223 if (it == lssRefUsrSel.end()) 00224 lssRefUsrSel.insert(refSelStrPair(grpID,refSelStrVec())); 00225 00226 lssRefUsrSel[grpID].push_back(refSelStr(segName, resID)); 00227 } 00228 00229 for (auto it=lssRefUsrSel.begin(); it!=lssRefUsrSel.end(); it++) { 00230 iout << iINFO << "LSS user selection for COM composition of group " 00231 << it->first << ":\n" << endi ; 00232 for (int i=0; i<it->second.size();i++) { 00233 iout << iINFO << "Segment " << it->second[i].segid 00234 << " ; residue " << it->second[i].resid << "\n" << endi ; 00235 } 00236 } 00237 } 00238 } 00239 00240 00241 00244 00245 00246 for (int atmInd = 0 ; atmInd < numAtoms; atmInd++) { 00247 00248 // If we are looking for QM-MM bonds, then we need to store extra info. 00249 if (simParams->qmBondOn) { 00250 00251 // We store both the qm group and the bond value 00252 if ( strcmp(simParams->qmColumn,"beta") == 0 ){ 00253 atmGrp = pdbP->atom(atmInd)->temperaturefactor() ; 00254 00255 bondVal = pdbP->atom(atmInd)->occupancy() ; 00256 } 00257 else { 00258 atmGrp = pdbP->atom(atmInd)->occupancy() ; 00259 00260 bondVal = pdbP->atom(atmInd)->temperaturefactor() ; 00261 } 00262 00263 qmBondValVec.push_back(bondVal); 00264 } 00265 else { 00266 00267 if ( strcmp(simParams->qmColumn,"beta") == 0 ){ 00268 atmGrp = pdbP->atom(atmInd)->temperaturefactor() ; 00269 } 00270 else { 00271 atmGrp = pdbP->atom(atmInd)->occupancy() ; 00272 } 00273 } 00274 00275 // We store all atom QM Group IDs in an array for 00276 // future transport to all PEs. 00277 qmAtomGroup[atmInd] = atmGrp; 00278 00279 // We store all atom's elements for quick access in the QM code. 00280 // Elements are used to tell the QM code the atom type. 00281 qmElementArray[atmInd] = new char[3]; 00282 strncpy(qmElementArray[atmInd],pdbP->atom(atmInd)->element(),3); 00283 00284 // For QM atoms, we keep global atom indices and parse different QM Group 00285 // IDs, keeping track of each groups size 00286 if (atmGrp > 0){ 00287 00288 if (simParams->fixedAtomsOn) { 00289 if ( fixedAtomFlags[atmInd] == 1 ) { 00290 iout << iERROR << "QM atom cannot be fixed in space!\n" 00291 << endi; 00292 NAMD_die("Error processing QM information."); 00293 } 00294 } 00295 00296 // Changes the VdW type of QM atoms. 00297 // This is may be used to counter balance the overpolarization 00298 // that QM atoms sufer. 00299 if (simParams->qmVDW) { 00300 // Creating a new type string 00301 std::string qmType(qmElementArray[atmInd]) ; 00302 // Erases empty spaces 00303 qmType.erase( 00304 std::remove_if(qmType.begin(), qmType.end(), isspace ), 00305 qmType.end()); 00306 // pre-pends a "q" to the element name. 00307 qmType = std::string("q") + qmType; 00308 00309 // iout << "QM VdW type: " << atoms[atmInd].vdw_type 00310 // << " atom type: " << atomNames[atmInd].atomtype 00311 // << " new type "<< qmType << "\n" << endi; 00312 00313 /* Look up the vdw constants for this atom */ 00314 // I am casting a non-const here because the function doesn't actually 00315 // change the string value, but it doesn't ask for a const char* as 00316 // an argument. 00317 params->assign_vdw_index(const_cast<char*>( qmType.c_str()), 00318 &(atoms[atmInd])); 00319 00320 // iout << "----> new VdW type: " << atoms[atmInd].vdw_type << "\n" << endi; 00321 } 00322 00323 // Indexes all global indices of QM atoms. 00324 qmAtmIndxVec.push_back(atmInd); 00325 00326 auto retVal = atmGrpSet.insert(atmGrp); 00327 00328 // If a new QM group is found, store the reverse reference in a map 00329 // and increase the total count. 00330 if (retVal.second) { 00331 // This mak makes the reverse identification from the group ID 00332 // to the sequential order in which each group was first found. 00333 qmGrpIDMap.insert(std::pair<BigReal,int>(atmGrp, qmNumGrps)); 00334 00335 // This vector keeps track of the group ID for each group 00336 qmGroupIDsVec.push_back(atmGrp); 00337 00338 // This counter is used to keep track of the sequential order in 00339 // which QM groups were first seen in the reference PDB file. 00340 qmNumGrps++ ; 00341 00342 // If this is a new QM group, initialize a new variable in the 00343 // vector to keep track of the number of atoms in this group. 00344 qmGrpSizeVec.push_back(1); 00345 00346 // For each new QM group, a new entry in the total charge 00347 // vector is created 00348 grpChrgVec.push_back( atoms[atmInd].charge ); 00349 } 00350 else { 00351 // If the QM group has already been seen, increase the count 00352 // of atoms in that group. 00353 qmGrpSizeVec[ qmGrpIDMap[atmGrp] ]++ ; 00354 00355 // If the QM group exists, adds current atom charge to total charge. 00356 grpChrgVec[ qmGrpIDMap[atmGrp] ] += atoms[atmInd].charge; 00357 } 00358 00359 // In case we are using live solvent selection 00360 if(simParams->qmLSSOn) { 00361 00362 int resID = pdbP->atom(atmInd)->residueseq(); 00363 char segName[5]; 00364 strncpy(segName, pdbP->atom(atmInd)->segmentname(),5); 00365 00366 int resSize = get_residue_size(segName,resID); 00367 00368 int i =-1, end =-1; 00369 00370 resLookup->lookup(segName,resID,&i, &end); 00371 00372 if ( strncasecmp(pdbP->atom(atmInd)->residuename(), 00373 simParams->qmLSSResname, 4) == 0) { 00374 00375 // We try to insert every residue from every atom 00376 qmSolvData *resP = lssGrpRes.find(qmSolvData(resID, i, resSize)); 00377 00378 if (resP != NULL) { 00379 resP->atmIDs.push_back(atmInd) ; 00380 // DebugM(3, "Existing residue " << resP->resID 00381 // << " from segID " << resP->segName 00382 // << " received atom " 00383 // << atmInd << "\n" ); 00384 } 00385 else { 00386 lssGrpRes.insert(qmSolvData(resID, i, resSize, segName, atmGrp)); 00387 // DebugM(3, lssGrpRes.size() << ") Inserted new residue: " 00388 // << resID << " from segID " << segName 00389 // << " with atom " 00390 // << i << "\n" ) ; 00391 } 00392 } 00393 else { 00394 // We store the QM atoms, per group, which are NOT part of 00395 // solvent molecules. 00396 00397 // Checks if we have a vector for this QM group. 00398 auto itNewGrp = lssGrpRefIDs.find(atmGrp) ; 00399 if (itNewGrp == lssGrpRefIDs.end()) { 00400 lssGrpRefIDs.insert(std::pair<Real, std::vector<int> >( 00401 atmGrp, std::vector<int>() ) ); 00402 } 00403 00404 switch (simParams->qmLSSMode) 00405 { 00406 00407 case QMLSSMODECOM: 00408 { 00409 // If we are using COM selection, checks if the atom 00410 // is part of the user selected residues 00411 for(int i=0; i<lssRefUsrSel[atmGrp].size(); i++) { 00412 if (lssRefUsrSel[atmGrp][i].resid == resID && 00413 strncmp(lssRefUsrSel[atmGrp][i].segid.c_str(), 00414 segName,5) == 0 ) { 00415 00416 lssGrpRefIDs[atmGrp].push_back(atmInd); 00417 totalNumRefAtms++; 00418 } 00419 } 00420 00421 } break; 00422 00423 case QMLSSMODEDIST: 00424 { 00425 00426 lssGrpRefIDs[atmGrp].push_back(atmInd); 00427 totalNumRefAtms++; 00428 00429 } break; 00430 00431 } 00432 } 00433 00434 } 00435 00436 } 00437 else if (atmGrp == 0) { 00438 00439 if(simParams->qmLSSOn) { 00440 00441 if ( strncasecmp(pdbP->atom(atmInd)->residuename(), 00442 simParams->qmLSSResname, 4) == 0) { 00443 00444 int resID = pdbP->atom(atmInd)->residueseq(); 00445 char segName[5]; 00446 strncpy(segName, pdbP->atom(atmInd)->segmentname(),5); 00447 00448 if (lssCurrClassResID < 0) { 00449 lssCurrClassResID = resID ; 00450 strncpy(lssCurrClassResSegID, segName,5); 00451 lssClassicResIndx = 0; 00452 } 00453 else if (lssCurrClassResID != resID || 00454 strcmp(lssCurrClassResSegID, segName) != 0 ) { 00455 lssCurrClassResID = resID ; 00456 strncpy(lssCurrClassResSegID, segName,5); 00457 lssClassicResIndx++; 00458 } 00459 00460 qmClassicSolv.insert(std::pair<int,int>(atmInd,lssClassicResIndx)); 00461 00462 } 00463 } 00464 } 00465 else if(atmGrp < 0) { 00466 iout << iERROR << "QM group ID cannot be less than zero!\n" << endi; 00467 NAMD_die("Error processing QM information."); 00468 } 00469 } 00470 00471 // Sanity check 00472 if (simParams->qmLSSOn) { 00473 if (lssGrpRes.size() == 0) 00474 NAMD_die("LSS was activated but there are no QM solvent molecules!\n"); 00475 00476 for (auto it=lssRefUsrSel.begin(); it!=lssRefUsrSel.end(); it++) { 00477 auto itTarget = qmGrpIDMap.find(it->first); 00478 if (itTarget == qmGrpIDMap.end()) { 00479 iout << iERROR << "QM group ID for LSS could not be found in input!" 00480 << " QM group ID: " << it->first << "\n" << endi; 00481 NAMD_die("Error processing QM information."); 00482 } 00483 } 00484 00485 DebugM(3,"We found " << lssClassicResIndx << " classical solvent molecules totalling " 00486 << qmClassicSolv.size() << " atoms.\n" ); 00487 00488 } 00489 00490 qmNumQMAtoms = qmAtmIndxVec.size(); 00491 00492 if (qmNumQMAtoms == 0) 00493 NAMD_die("No QM atoms were found in this QM simulation!") ; 00494 00495 iout << iINFO << "Number of QM atoms (excluding Dummy atoms): " << 00496 qmNumQMAtoms << "\n" << endi; 00497 00498 qmAtmChrg = new Real[qmNumQMAtoms] ; 00499 qmAtmIndx = new int[qmNumQMAtoms] ; 00500 for (int i=0; i<qmNumQMAtoms; i++) { 00501 // qmAtmChrg gets initialized with the PSF charges at the end of this 00502 // function, but values may change as QM steps are taken. 00503 qmAtmChrg[i] = 0; 00504 qmAtmIndx[i] = qmAtmIndxVec[i] ; 00505 } 00506 00507 // This map relates the QM group index with a vector of pairs 00508 // of bonded MM-QM atoms (with the bonded atoms ins this order: 00509 // MM first, QM second). 00510 std::map<int, std::vector<std::pair<int,int> > > qmGrpIDBonds ; 00511 int bondCounter = 0; 00512 if (simParams->qmBondOn) { 00513 00514 // Checks all bonds for QM-MM pairs. 00515 // Makes sanity checks against QM-QM pairs and MM-MM pairs that 00516 // are flagged by the user to be bonds between QM and MM regions. 00517 // QM-QM bonds will be removed in another function. 00518 for (int bndIt = 0; bndIt < numBonds; bndIt++) { 00519 00520 bond curr = bonds[bndIt] ; 00521 00522 // In case either atom has a non-zero 00523 if ( qmBondValVec[curr.atom1] != 0 && 00524 qmBondValVec[curr.atom2] != 0 ) { 00525 00526 // Sanity checks (1 of 2) 00527 if (qmAtomGroup[curr.atom1] != 0 && 00528 qmAtomGroup[curr.atom2] != 0) { 00529 iout << iERROR << "Atoms " << curr.atom1 << " and " << 00530 curr.atom2 << " are assigned as QM atoms.\n" << endi; 00531 NAMD_die("Error in QM-MM bond assignment.") ; 00532 } 00533 00534 // Sanity checks (2 of 2) 00535 if (qmAtomGroup[curr.atom1] == 0 && 00536 qmAtomGroup[curr.atom2] == 0) { 00537 iout << iERROR << "Atoms " << curr.atom1 << " and " << 00538 curr.atom2 << " are assigned as MM atoms.\n" << endi; 00539 NAMD_die("Error in QM-MM bond assignment.") ; 00540 } 00541 00542 int currGrpID = 0; 00543 std::pair<int,int> newPair(0,0); 00544 00545 // We create a QM-MM pair with the MM atom first 00546 if (qmAtomGroup[curr.atom1] != 0) { 00547 newPair = std::pair<int,int>(curr.atom2, curr.atom1) ; 00548 currGrpID = qmAtomGroup[curr.atom1]; 00549 } else { 00550 newPair = std::pair<int,int>(curr.atom1, curr.atom2) ; 00551 currGrpID = qmAtomGroup[curr.atom2]; 00552 } 00553 00554 int grpIndx = qmGrpIDMap[currGrpID] ; 00555 00556 // Stores bonds in vectors key-ed by the QM group ID of the QM atom. 00557 auto retIter = qmGrpIDBonds.find(grpIndx) ; 00558 // In case thi QM-MM bonds belong to a QM group we have not seen 00559 // yet, stores a new vector in the map. 00560 if (retIter == qmGrpIDBonds.end()) { 00561 qmGrpIDBonds.insert(std::pair<BigReal,std::vector<std::pair<int,int> > >( 00562 grpIndx, std::vector<std::pair<int,int> >() ) ); 00563 } 00564 00565 qmGrpIDBonds[grpIndx].push_back( newPair ); 00566 00567 bondCounter++ ; 00568 } 00569 00570 00571 } 00572 00573 // iout << "Finished processing QM-MM bonds.\n" << endi ; 00574 00575 if (bondCounter == 0) 00576 iout << iWARN << "We found " << bondCounter << " QM-MM bonds.\n" << endi ; 00577 else 00578 iout << iINFO << "We found " << bondCounter << " QM-MM bonds.\n" << endi ; 00579 00580 } 00581 00582 // Initializes several arrays used to setup the QM simulation. 00583 qmNumBonds = bondCounter ; 00584 00585 qmMMBond = new int*[qmNumBonds]; 00586 00587 qmDummyBondVal = new BigReal[qmNumBonds]; 00588 00589 qmMMChargeTarget = new int*[qmNumBonds] ; 00590 qmMMNumTargs = new int[qmNumBonds] ; 00591 00592 qmDummyElement = new char*[qmNumBonds] ; 00593 00594 00595 qmNumGrps = atmGrpSet.size(); 00596 00597 qmGrpSizes = new int[qmNumGrps]; 00598 00599 qmGrpID = new Real[qmNumGrps]; 00600 00601 qmGrpChrg = new Real[qmNumGrps]; 00602 00603 qmGrpMult = new Real[qmNumGrps] ; 00604 00605 qmGrpNumBonds = new int[qmNumGrps]; 00606 00607 qmGrpBonds = new int*[qmNumGrps]; 00608 qmMMBondedIndx = new int*[qmNumGrps]; 00609 00610 00613 00614 00615 // We first initialize the multiplicity vector with 00616 // default values, then populate it with user defined values. 00617 for (int grpIter = 0; grpIter < qmNumGrps; grpIter++) { 00618 qmGrpMult[grpIter] = 0; 00619 } 00620 00621 int multCount = 0; 00622 StringList *current = cfgList->find("QMMult"); 00623 for ( ; current; current = current->next ) { 00624 00625 auto strVec = split( std::string(current->data) , " "); 00626 00627 if (strVec.size() != 2 ) { 00628 iout << iERROR << "Format error in QM Multiplicity string: " 00629 << current->data 00630 << "\n" << endi; 00631 NAMD_die("Error processing QM information."); 00632 } 00633 00634 std::stringstream storConv ; 00635 00636 storConv << strVec[0] ; 00637 Real grpID ; 00638 storConv >> grpID; 00639 if (storConv.fail()) { 00640 iout << iERROR << "Error parsing QMMult selection: " 00641 << current->data 00642 << "\n" << endi; 00643 NAMD_die("Error processing QM information."); 00644 } 00645 00646 storConv.clear() ; 00647 storConv << strVec[1]; 00648 Real multiplicity ; 00649 storConv >> multiplicity; 00650 if (storConv.fail()) { 00651 iout << iERROR << "Error parsing QMMult selection: " 00652 << current->data 00653 << "\n" << endi; 00654 NAMD_die("Error processing QM information."); 00655 } 00656 00657 auto it = qmGrpIDMap.find(grpID); 00658 00659 if (it == qmGrpIDMap.end()) { 00660 iout << iERROR << "Error parsing QMMult selection. Could not find QM group ID: " 00661 << grpID 00662 << "\n" << endi; 00663 NAMD_die("Error processing QM information."); 00664 } 00665 else { 00666 iout << iINFO << "Applying user defined multiplicity " 00667 << multiplicity << " to QM group ID " << grpID << "\n" << endi; 00668 qmGrpMult[it->second] = multiplicity; 00669 } 00670 00671 multCount++; 00672 } 00673 00674 if (multCount != qmNumGrps && simParams->qmFormat == QMFormatORCA) { 00675 NAMD_die("ORCA neds multiplicity values for all QM regions."); 00676 } 00677 00678 00681 00682 00683 for (size_t grpIndx = 0; grpIndx < qmGrpSizeVec.size(); grpIndx++) { 00684 00685 bool nonInteger = true; 00686 if ((fabsf(roundf(grpChrgVec[grpIndx]) - grpChrgVec[grpIndx]) <= 0.001f) ) { 00687 grpChrgVec[grpIndx] = roundf(grpChrgVec[grpIndx]) ; 00688 nonInteger = false; 00689 } 00690 00691 iout << iINFO << grpIndx + 1 << ") Group ID: " << qmGroupIDsVec[grpIndx] 00692 << " ; Group size: " << qmGrpSizeVec[grpIndx] << " atoms" 00693 << " ; Total PSF charge: " << grpChrgVec[grpIndx] << "\n" << endi ; 00694 00695 if (nonInteger && simParams->PMEOn) 00696 NAMD_die("QM atoms do not add up to a whole charge, which is needed for PME.") ; 00697 } 00698 00699 int chrgCount = 0; 00700 current = cfgList->find("QMCharge"); 00701 for ( ; current; current = current->next ) { 00702 00703 auto strVec = split( std::string(current->data) , " "); 00704 00705 if (strVec.size() != 2 ) { 00706 iout << iERROR << "Format error in QM Charge string: " 00707 << current->data 00708 << "\n" << endi; 00709 NAMD_die("Error processing QM information."); 00710 } 00711 00712 std::stringstream storConv ; 00713 00714 storConv << strVec[0] ; 00715 Real grpID ; 00716 storConv >> grpID; 00717 if (storConv.fail()) { 00718 iout << iERROR << "Error parsing QMCharge selection: " 00719 << current->data 00720 << "\n" << endi; 00721 NAMD_die("Error processing QM information."); 00722 } 00723 00724 storConv.clear() ; 00725 storConv << strVec[1]; 00726 Real charge ; 00727 storConv >> charge; 00728 if (storConv.fail()) { 00729 iout << iERROR << "Error parsing QMCharge selection: " 00730 << current->data 00731 << "\n" << endi; 00732 NAMD_die("Error processing QM information."); 00733 } 00734 00735 auto it = qmGrpIDMap.find(grpID); 00736 00737 if (it == qmGrpIDMap.end()) { 00738 iout << iERROR << "Error parsing QMCharge selection. Could not find QM group ID: " 00739 << grpID 00740 << "\n" << endi; 00741 NAMD_die("Error processing QM information."); 00742 } 00743 else { 00744 iout << iINFO << "Found user defined charge " 00745 << charge << " for QM group ID " << grpID << ". Will ignore PSF charge.\n" << endi; 00746 grpChrgVec[it->second] = charge; 00747 } 00748 00749 chrgCount++; 00750 } 00751 00752 simParams->qmMOPACAddConfigChrg = false; 00753 // Checks if QM group charges can be defined for MOPAC. 00754 // Since charge can be defined in QMConfigLine, we need extra logic. 00755 if (simParams->qmFormat == QMFormatMOPAC) { 00756 00757 // Checks if group charge was supplied in config line for MOPAC. 00758 std::string::size_type chrgLoc = std::string::npos ; 00759 00760 // This will hold a sting with the first user supplied configuration line for MOPAC. 00761 std::string configLine(cfgList->find("QMConfigLine")->data) ; 00762 chrgLoc = configLine.find("CHARGE") ; 00763 00764 if ( chrgLoc != std::string::npos ) { 00765 iout << iINFO << "Found user defined charge in command line. This \ 00766 will be used for all QM regions and will take precedence over all other charge \ 00767 definitions.\n" << endi; 00768 } 00769 else if ( chrgLoc == std::string::npos && (simParams->qmChrgFromPSF || chrgCount == qmNumGrps) ) { 00770 // If no charge was defined in the configuration line, gets from PSF and/or 00771 // from user defined charges (through the QMCharge keyword). 00772 simParams->qmMOPACAddConfigChrg = true; 00773 } 00774 else 00775 { 00776 // If we could nont find a charge definition in the config line AND 00777 // no specific charge was selected for each QM region through QMCharge AND 00778 // the QMChargeFromPSF was not turned ON, the we stop NAMD and scream at the user. 00779 // if ( chrgLoc == std::string::npos && (chrgCount != qmNumGrps ) && !simParams->qmChrgFromPSF) 00780 iout << iERROR << "Could not find charge for all QM groups. For MOPAC, \ 00781 charges can be defined through QMConfigLine, QMCharge or QMChargeFromPSF keywords.\n" << endi; 00782 NAMD_die("Error processing QM information."); 00783 } 00784 00785 } 00786 00787 if (simParams->qmFormat == QMFormatORCA ) { 00788 if ((chrgCount != qmNumGrps ) && !simParams->qmChrgFromPSF) { 00789 // If we are not supposed to get charges from the PSF, and not 00790 // enough charges were set to cover all QM regions, 00791 // we stop NAMD and scream at the user. 00792 iout << iERROR << "Could not find charge for all QM groups. For ORCA, \ 00793 charges can be defined through QMCharge or QMChargeFromPSF keywords.\n" << endi; 00794 NAMD_die("Error processing QM information."); 00795 } 00796 } 00797 00798 // If mechanichal embedding was requested but we have QM-MM bonds, we need 00799 // to send extra info to ComputeQM to preserve calculation speed. 00800 if (qmNumBonds > 0 && qmNoPC) { 00801 qmMeNumBonds = qmNumBonds; 00802 qmMeMMindx = new int[qmMeNumBonds] ; 00803 qmMeQMGrp = new Real[qmMeNumBonds] ; 00804 } 00805 else { 00806 qmMeNumBonds = 0 ; 00807 qmMeMMindx = NULL; 00808 qmMeQMGrp = NULL; 00809 } 00810 00811 00814 00815 00816 bondCounter = 0; 00817 for (int grpIter = 0; grpIter < qmNumGrps; grpIter++) { 00818 00819 qmGrpID[grpIter] = qmGroupIDsVec[grpIter] ; 00820 00821 qmGrpSizes[grpIter] = qmGrpSizeVec[grpIter] ; 00822 00823 qmGrpChrg[grpIter] = grpChrgVec[grpIter]; 00824 00825 // iout << "Loaded " << qmGrpSizes[grpIter] << " atoms to this group.\n" << endi; 00826 00827 int currNumbBonds = qmGrpIDBonds[grpIter].size() ; 00828 00829 // Assigns the number of bonds that the current QM group has. 00830 qmGrpNumBonds[grpIter] = currNumbBonds; 00831 00832 if (currNumbBonds > 0) { 00833 00834 qmGrpBonds[grpIter] = new int[currNumbBonds]; 00835 qmMMBondedIndx[grpIter] = new int[currNumbBonds]; 00836 00837 for (int bndIter = 0; bndIter<currNumbBonds; bndIter++) { 00838 00839 // Adds the bonds to the overall sequential list. 00840 qmMMBond[bondCounter] = new int[2] ; 00841 qmMMBond[bondCounter][0] = qmGrpIDBonds[grpIter][bndIter].first ; 00842 qmMMBond[bondCounter][1] = qmGrpIDBonds[grpIter][bndIter].second ; 00843 00844 // For the current QM group, and the current bond, gets the bond index. 00845 qmGrpBonds[grpIter][bndIter] = bondCounter; 00846 00847 // For the current QM group, and the current bond, gets the MM atom. 00848 qmMMBondedIndx[grpIter][bndIter] = qmGrpIDBonds[grpIter][bndIter].first ; 00849 00850 // Assign the default value of dummy element 00851 qmDummyElement[bondCounter] = new char[3]; 00852 strcpy(qmDummyElement[bondCounter],"H\0"); 00853 00854 // Determines the distance that will separate the new Dummy atom 00855 // and the Qm atom to which it will be bound. 00856 bondVal = 0; 00857 if (simParams->qmBondDist) { 00858 if (qmBondValVec[qmMMBond[bondCounter][0]] != 00859 qmBondValVec[qmMMBond[bondCounter][1]] 00860 ) { 00861 iout << iERROR << "qmBondDist is ON but the values in the bond column are different for atom " 00862 << qmMMBond[bondCounter][0] << " and " << qmMMBond[bondCounter][1] 00863 << ".\n" << endi ; 00864 NAMD_die("Error in QM-MM bond processing."); 00865 } 00866 00867 bondVal = qmBondValVec[qmMMBond[bondCounter][0]] ; 00868 } 00869 else { 00870 00871 if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"C") == 0 ) { 00872 bondVal = 1.09 ; 00873 } 00874 else if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"O") == 0 ) { 00875 bondVal = 0.98 ; 00876 } 00877 else if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"N") == 0 ) { 00878 bondVal = 0.99 ; 00879 } 00880 else { 00881 iout << iERROR << "qmBondDist is OFF but the QM atom has no default distance value. Atom " 00882 << qmMMBond[bondCounter][1] << ", with element: " 00883 << qmElementArray[qmMMBond[bondCounter][1]] 00884 << ".\n" << endi ; 00885 NAMD_die("Error in QM-MM bond processing."); 00886 } 00887 00888 } 00889 00890 iout << iINFO << "MM-QM pair: " << qmMMBond[bondCounter][0] << ":" 00891 << qmMMBond[bondCounter][1] 00892 << " -> Value (distance or ratio): " << bondVal 00893 << " (QM Group " << grpIter << " ID " << qmGrpID[grpIter] << ")" 00894 << "\n" << endi ; 00895 00896 qmDummyBondVal[bondCounter] = bondVal; 00897 00898 // In case we are preparing for a mechanical embedding simulation 00899 // with no point charges, populate the following vectors 00900 if (qmMeNumBonds > 0) { 00901 qmMeMMindx[bondCounter] = qmMMBond[bondCounter][0]; 00902 qmMeQMGrp[bondCounter] = qmGrpID[grpIter]; 00903 } 00904 00905 bondCounter++ ; 00906 } 00907 } 00908 00909 } 00910 00911 00914 00915 current = NULL; 00916 if (qmNumBonds > 0) 00917 current = cfgList->find("QMLinkElement"); 00918 00919 int numParsedLinkElem = 0; 00920 for ( ; current != NULL; current = current->next ) { 00921 00922 DebugM(3,"Parsing link atom element: " << current->data << "\n" ); 00923 00924 strVec = split( std::string(current->data) , " "); 00925 00926 // We need the two atoms that compose the QM-MM bonds and 00927 // then the element. 00928 if (strVec.size() != 3 && qmNumBonds > 1) { 00929 iout << iERROR << "Format error in QM link atom element selection: " 00930 << current->data 00931 << "\n" << endi; 00932 NAMD_die("Error processing QM information."); 00933 } 00934 00935 // If there is only one QM-MM bond, we can accept the element only. 00936 if (strVec.size() != 1 && strVec.size() != 3 && qmNumBonds == 1) { 00937 iout << iERROR << "Format error in QM link atom element selection: " 00938 << current->data 00939 << "\n" << endi; 00940 NAMD_die("Error processing QM information."); 00941 } 00942 00943 std::stringstream storConv ; 00944 int bondAtom1, bondAtom2; 00945 std::string linkElement ; 00946 00947 if (strVec.size() == 1) { 00948 linkElement = strVec[0].substr(0,2); 00949 } 00950 else if (strVec.size() == 3) { 00951 00952 storConv << strVec[0] ; 00953 storConv >> bondAtom1; 00954 if (storConv.fail()) { 00955 iout << iERROR << "Error parsing link atom element selection: " 00956 << current->data 00957 << "\n" << endi; 00958 NAMD_die("Error processing QM information."); 00959 } 00960 00961 storConv.clear() ; 00962 storConv << strVec[1]; 00963 storConv >> bondAtom2; 00964 if (storConv.fail()) { 00965 iout << iERROR << "Error parsing link atom element selection: " 00966 << current->data 00967 << "\n" << endi; 00968 NAMD_die("Error processing QM information."); 00969 } 00970 00971 linkElement = strVec[2].substr(0,2); 00972 } 00973 00974 numParsedLinkElem++; 00975 00976 if (numParsedLinkElem>qmNumBonds) { 00977 iout << iERROR << "More link atom elements were set than there" 00978 " are QM-MM bonds.\n" << endi; 00979 NAMD_die("Error processing QM information."); 00980 } 00981 00982 int bondIter; 00983 00984 if (strVec.size() == 1) { 00985 bondIter = 0; 00986 } 00987 else if (strVec.size() == 3) { 00988 00989 Bool foundBond = false; 00990 00991 for (bondIter=0; bondIter<qmNumBonds; bondIter++) { 00992 00993 if ( (qmMMBond[bondIter][0] == bondAtom1 && 00994 qmMMBond[bondIter][1] == bondAtom2 ) || 00995 (qmMMBond[bondIter][0] == bondAtom2 && 00996 qmMMBond[bondIter][1] == bondAtom1 ) ) { 00997 00998 foundBond = true; 00999 break; 01000 } 01001 } 01002 01003 if (! foundBond) { 01004 iout << iERROR << "Error parsing link atom element selection: " 01005 << current->data 01006 << "\n" << endi; 01007 iout << iERROR << "Atom pair was not found to be a QM-MM bond.\n" 01008 << endi; 01009 NAMD_die("Error processing QM information."); 01010 } 01011 } 01012 01013 strcpy(qmDummyElement[bondIter],linkElement.c_str()); 01014 qmDummyElement[bondIter][2] = '\0'; 01015 01016 iout << iINFO << "Applying user defined link atom element " 01017 << qmDummyElement[bondIter] << " to QM-MM bond " 01018 << bondIter << ": " << qmMMBond[bondIter][0] 01019 << " - " << qmMMBond[bondIter][1] 01020 << "\n" << endi; 01021 } 01022 01023 01024 01027 01028 01029 int32 **bondsWithAtomLocal = NULL ; 01030 int32 *byAtomSizeLocal = NULL; 01031 ObjectArena <int32 >* tmpArenaLocal = NULL; 01032 if (qmNumBonds > 0) { 01033 01034 bondsWithAtomLocal = new int32 *[numAtoms]; 01035 byAtomSizeLocal = new int32[numAtoms]; 01036 tmpArenaLocal = new ObjectArena<int32>; 01037 01038 // Build the bond lists 01039 for (int i=0; i<numAtoms; i++) 01040 { 01041 byAtomSizeLocal[i] = 0; 01042 } 01043 for (int i=0; i<numRealBonds; i++) 01044 { 01045 byAtomSizeLocal[bonds[i].atom1]++; 01046 byAtomSizeLocal[bonds[i].atom2]++; 01047 } 01048 for (int i=0; i<numAtoms; i++) 01049 { 01050 bondsWithAtomLocal[i] = tmpArenaLocal->getNewArray(byAtomSizeLocal[i]+1); 01051 bondsWithAtomLocal[i][byAtomSizeLocal[i]] = -1; 01052 byAtomSizeLocal[i] = 0; 01053 } 01054 for (int i=0; i<numRealBonds; i++) 01055 { 01056 int a1 = bonds[i].atom1; 01057 int a2 = bonds[i].atom2; 01058 bondsWithAtomLocal[a1][byAtomSizeLocal[a1]++] = i; 01059 bondsWithAtomLocal[a2][byAtomSizeLocal[a2]++] = i; 01060 } 01061 } 01062 01063 // In this loops we try to find other bonds in which the MM atoms from 01064 // QM-MM bonds may be involved. The other MM atoms (which we will call M2 and M3) 01065 // will be involved in charge manipulation. See ComputeQM.C for details. 01066 for (int qmBndIt = 0; qmBndIt < qmNumBonds; qmBndIt++) { 01067 01068 // The charge targets are accumulated in a temporary vector and then 01069 // transfered to an array that will be transmited to the ComputeQMMgr object. 01070 std::vector<int> chrgTrgt ; 01071 01072 int MM1 = qmMMBond[qmBndIt][0], MM2, MM2BondIndx, MM3, MM3BondIndx; 01073 01074 switch (simParams->qmBondScheme) { 01075 01076 case QMSCHEMERCD: 01077 01078 case QMSCHEMECS: 01079 { 01080 // Selects ALL MM2 atoms. 01081 for (int i=0; i<byAtomSizeLocal[MM1]; i++) { 01082 MM2BondIndx = bondsWithAtomLocal[MM1][i] ; 01083 01084 // Checks which of the atoms in the bond structure is the 01085 // MM2 atom. 01086 if (bonds[MM2BondIndx].atom1 == MM1) 01087 MM2 = bonds[MM2BondIndx].atom2; 01088 else 01089 MM2 = bonds[MM2BondIndx].atom1; 01090 01091 // In case the bonded atom is a QM atom, 01092 // skips the index. 01093 if (qmAtomGroup[MM2] > 0) 01094 continue; 01095 01096 chrgTrgt.push_back(MM2); 01097 } 01098 01099 } break; 01100 01101 case QMSCHEMEZ3: 01102 { 01103 // Selects all MM3 atoms 01104 for (int i=0; i<byAtomSizeLocal[MM1]; i++) { 01105 MM2BondIndx = bondsWithAtomLocal[MM1][i] ; 01106 01107 // Checks which of the atoms in the bond structure is the 01108 // MM2 atom. 01109 if (bonds[MM2BondIndx].atom1 == MM1) 01110 MM2 = bonds[MM2BondIndx].atom2; 01111 else 01112 MM2 = bonds[MM2BondIndx].atom1; 01113 01114 // In case the bonded atom is a QM atom, 01115 // skips the index. 01116 if (qmAtomGroup[MM2] > 0) 01117 continue; 01118 01119 for (int i=0; i<byAtomSizeLocal[MM2]; i++) { 01120 MM3BondIndx = bondsWithAtomLocal[MM2][i] ; 01121 01122 // Checks which of the atoms in the bond structure is the 01123 // MM3 atom. 01124 if (bonds[MM3BondIndx].atom1 == MM2) 01125 MM3 = bonds[MM3BondIndx].atom2; 01126 else 01127 MM3 = bonds[MM3BondIndx].atom1; 01128 01129 // In case the bonded atom is a QM atom, 01130 // skips the index. 01131 // We also keep the search from going back to MM1. 01132 if (qmAtomGroup[MM3] > 0 || MM3 == MM1) 01133 continue; 01134 01135 chrgTrgt.push_back(MM3); 01136 } 01137 01138 } 01139 01140 }; 01141 01142 case QMSCHEMEZ2: 01143 { 01144 // Selects all MM2 atoms 01145 for (int i=0; i<byAtomSizeLocal[MM1]; i++) { 01146 MM2BondIndx = bondsWithAtomLocal[MM1][i] ; 01147 01148 // Checks which of the atoms in the bond structure is the 01149 // MM2 atom. 01150 if (bonds[MM2BondIndx].atom1 == MM1) 01151 MM2 = bonds[MM2BondIndx].atom2; 01152 else 01153 MM2 = bonds[MM2BondIndx].atom1; 01154 01155 // In case the bonded atom is a QM atom, 01156 // skips the index. 01157 if (qmAtomGroup[MM2] > 0) 01158 continue; 01159 01160 chrgTrgt.push_back(MM2); 01161 } 01162 01163 }; 01164 01165 case QMSCHEMEZ1: 01166 { 01167 // Selects all MM1 atoms 01168 chrgTrgt.push_back(MM1); 01169 } break; 01170 } 01171 01172 01173 qmMMChargeTarget[qmBndIt] = new int[chrgTrgt.size()] ; 01174 qmMMNumTargs[qmBndIt] = chrgTrgt.size(); 01175 01176 DebugM(3, "MM-QM bond " << qmBndIt << "; MM atom " 01177 << qmMMBond[qmBndIt][0] << " conections: \n" ); 01178 01179 for (size_t i=0; i < chrgTrgt.size(); i++) { 01180 qmMMChargeTarget[qmBndIt][i] = chrgTrgt[i]; 01181 DebugM(3,"MM Bonded to: " << chrgTrgt[i] << "\n" ); 01182 } 01183 01184 chrgTrgt.clear(); 01185 } 01186 01187 if (bondsWithAtomLocal != NULL) 01188 delete [] bondsWithAtomLocal; bondsWithAtomLocal = 0; 01189 if (byAtomSizeLocal != NULL) 01190 delete [] byAtomSizeLocal; byAtomSizeLocal = 0; 01191 if (tmpArenaLocal != NULL) 01192 delete tmpArenaLocal; tmpArenaLocal = 0; 01193 01194 01197 01198 01199 if(simParams->qmLSSOn) { 01200 01201 std::map<Real,int> grpLSSSize ; 01202 std::map<Real,int>::iterator itGrpSize; 01203 01204 qmLSSTotalNumAtms = 0; 01205 qmLSSResidueSize = 0; 01206 01207 if (simParams->qmLSSFreq == 0) 01208 qmLSSFreq = simParams->stepsPerCycle ; 01209 else 01210 qmLSSFreq = simParams->qmLSSFreq; 01211 01212 #ifdef DEBUG_QM 01213 int resSize = -1; 01214 #endif 01215 01216 std::map<Real, int> grpNumLSSRes; 01217 std::map<Real, int>::iterator itGrpNumRes; 01218 01219 for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) { 01220 01221 if (it->atmIDs.size() != it->size) { 01222 iout << iERROR << "The number of atoms loaded for residue " 01223 << it->resID << " does not match the expected for this residue type.\n" 01224 << endi; 01225 NAMD_die("Error parsing data for LSS."); 01226 } 01227 01228 qmLSSTotalNumAtms += it->size; 01229 01230 #ifdef DEBUG_QM 01231 if (resSize < 0) resSize = it->size ; 01232 if (resSize > 0 and resSize != it->size) { 01233 iout << iERROR << "The number of atoms loaded for residue " 01234 << it->resID << " does not match previously loaded residues.\n" 01235 << endi; 01236 NAMD_die("Error parsing data for LSS."); 01237 } 01238 01239 // DebugM(3,"Residue " << it->resID << ": " << it->segName 01240 // << " - from " << it->begAtmID << " with size " 01241 // << it->size << " (QM ID: " << it->qmGrpID 01242 // << ") has " << it->atmIDs.size() << " atoms: \n" ) ; 01243 // for (int i=0; i<it->atmIDs.size(); i++) 01244 // DebugM(3, it->atmIDs[i] << "\n" ); 01245 #endif 01246 01247 // Calculating total number of atoms per group 01248 itGrpSize = grpLSSSize.find(it->qmGrpID) ; 01249 if (itGrpSize != grpLSSSize.end()) 01250 itGrpSize->second += it->size; 01251 else 01252 grpLSSSize.insert(std::pair<Real,int>(it->qmGrpID, it->size)); 01253 01254 // Calculating total number of solvent residues per group 01255 itGrpNumRes = grpNumLSSRes.find(it->qmGrpID) ; 01256 if (itGrpNumRes != grpNumLSSRes.end()) 01257 itGrpNumRes->second += 1; 01258 else 01259 grpNumLSSRes.insert(std::pair<Real,int>(it->qmGrpID, 1)); 01260 } 01261 01262 qmLSSResidueSize = lssGrpRes.begin()->size ; 01263 01264 qmLSSSize = new int[qmNumGrps]; 01265 01266 qmLSSIdxs = new int[qmLSSTotalNumAtms]; 01267 int lssAtmIndx = 0; 01268 01269 switch (simParams->qmLSSMode) { 01270 01271 case QMLSSMODECOM: 01272 { 01273 01274 qmLSSRefSize = new int[qmNumGrps]; 01275 01276 qmLSSMass = new Mass[qmLSSTotalNumAtms]; 01277 01278 qmLSSRefIDs = new int[totalNumRefAtms]; 01279 qmLSSRefMass = new Mass[totalNumRefAtms]; 01280 int lssRefIndx = 0; 01281 01282 for (int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) { 01283 01284 itGrpSize = grpNumLSSRes.find(qmGrpID[grpIndx]) ; 01285 01286 if (itGrpSize != grpNumLSSRes.end()) 01287 qmLSSSize[grpIndx] = itGrpSize->second; 01288 else 01289 qmLSSSize[grpIndx] = 0; 01290 01291 for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) { 01292 01293 if (it->qmGrpID == qmGrpID[grpIndx]) { 01294 for (int i=0; i<it->atmIDs.size(); i++) { 01295 qmLSSIdxs[lssAtmIndx] = it->atmIDs[i]; 01296 qmLSSMass[lssAtmIndx] = atoms[it->atmIDs[i]].mass; 01297 lssAtmIndx++; 01298 } 01299 } 01300 } 01301 01302 DebugM(4, "QM group " << qmGrpID[grpIndx] 01303 << " has " << qmLSSSize[grpIndx] << " solvent molecules, " 01304 << "totalling " << grpLSSSize[qmGrpID[grpIndx]] 01305 << " atoms (check " << lssAtmIndx << ").\n" ); 01306 01307 qmLSSRefSize[grpIndx] = lssGrpRefIDs[qmGrpID[grpIndx]].size(); 01308 for(int i=0; i < qmLSSRefSize[grpIndx]; i++) { 01309 qmLSSRefIDs[lssRefIndx] = lssGrpRefIDs[qmGrpID[grpIndx]][i]; 01310 qmLSSRefMass[lssRefIndx] = atoms[qmLSSRefIDs[lssRefIndx]].mass; 01311 lssRefIndx++; 01312 } 01313 01314 DebugM(4, "QM group " << qmGrpID[grpIndx] 01315 << " has " << qmLSSRefSize[grpIndx] << " non-solvent atoms for LSS.\n" ); 01316 } 01317 01318 } break ; 01319 01320 case QMLSSMODEDIST: 01321 { 01322 for (int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) { 01323 01324 itGrpSize = grpNumLSSRes.find(qmGrpID[grpIndx]) ; 01325 01326 if (itGrpSize != grpNumLSSRes.end()) 01327 qmLSSSize[grpIndx] = itGrpSize->second; 01328 else 01329 qmLSSSize[grpIndx] = 0; 01330 01331 for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) { 01332 01333 if (it->qmGrpID == qmGrpID[grpIndx]) { 01334 for (int i=0; i<it->atmIDs.size(); i++) { 01335 qmLSSIdxs[lssAtmIndx] = it->atmIDs[i]; 01336 lssAtmIndx++; 01337 } 01338 } 01339 } 01340 01341 DebugM(4, "QM group " << qmGrpID[grpIndx] 01342 << " has " << qmLSSSize[grpIndx] << " solvent molecules, " 01343 << "totalling " << grpLSSSize[qmGrpID[grpIndx]] 01344 << " atoms (check " << lssAtmIndx << ").\n" ); 01345 } 01346 01347 } break ; 01348 01349 } 01350 } 01351 01352 01355 01356 01357 PDB *customPCPDB; 01358 01359 // In case we have a custom and fixed set of point charges for each QM group, 01360 // we process the files containing information. 01361 current = NULL; 01362 if (simParams->qmCustomPCSel) { 01363 current = cfgList->find("QMCustomPCFile"); 01364 } 01365 01366 std::map<Real,std::vector<int> > qmPCVecMap ; 01367 01368 int numParsedPBDs = 0; 01369 for ( ; current != NULL; current = current->next ) { 01370 01371 iout << iINFO << "Parsing QM Custom PC file " << current->data << "\n" << endi; 01372 customPCPDB = new PDB(current->data); 01373 01374 if (customPCPDB->num_atoms() != numAtoms) 01375 NAMD_die("Number of atoms in QM Custom PC PDB file doesn't match coordinate PDB"); 01376 01377 std::vector< int > currPCSel ; 01378 Real currQMID = 0 ; 01379 int currGrpSize = 0 ; 01380 01381 for (int atmInd = 0 ; atmInd < numAtoms; atmInd++) { 01382 01383 BigReal beta = customPCPDB->atom(atmInd)->temperaturefactor() ; 01384 BigReal occ = customPCPDB->atom(atmInd)->occupancy() ; 01385 01386 if ( beta != 0 && occ != 0) 01387 NAMD_die("An atom cannot be marked as QM and poitn charge simultaneously!"); 01388 01389 // If this is not a QM atom and 01390 if (occ != 0) { 01391 currPCSel.push_back(atmInd) ; 01392 } 01393 01394 if (beta != 0) { 01395 if (pdbP->atom(atmInd)->temperaturefactor() != beta) 01396 NAMD_die("QM Group selection is different in reference PDB and Custom-PC PDB!"); 01397 01398 if (currQMID == 0) { 01399 // If this is the first QM atom we find, keep the QM Group ID. 01400 currQMID = beta; 01401 } 01402 else { 01403 // For all other atoms, check if it is the same group. It must be!! 01404 if (currQMID != beta) 01405 NAMD_die("Found two different QM group IDs in this file!"); 01406 } 01407 01408 currGrpSize++; 01409 } 01410 01411 } 01412 01413 if (currGrpSize != qmGrpSizeVec[ qmGrpIDMap[currQMID] ]) 01414 NAMD_die("Reference PDB and Custom-PC PDB have different QM group sizes!") ; 01415 01416 qmPCVecMap.insert(std::pair<Real,std::vector<int> >( 01417 currQMID, currPCSel )); 01418 01419 numParsedPBDs++; 01420 delete customPCPDB; 01421 } 01422 01423 delete pdbP; 01424 01425 if (numParsedPBDs != qmNumGrps && simParams->qmCustomPCSel) { 01426 iout << iWARN << "The number of files provided for custom point " 01427 "charges is not equal to the number of QM groups!\n" << endi; 01428 } 01429 01430 // Initializes an array with the number of Custom Point Charges per 01431 // QM group. 01432 qmCustPCSizes = new int[qmNumGrps]; 01433 for (int i=0; i<qmNumGrps; i++) 01434 qmCustPCSizes[i] = 0; 01435 01436 qmTotCustPCs = 0; 01437 01438 // Stores the size of each Custom PC vector in the array. 01439 // We may not have PCs for all QM groups. 01440 for (auto it = qmPCVecMap.begin(); it != qmPCVecMap.end(); it++) { 01441 qmTotCustPCs += (*it).second.size(); 01442 int qmIndx = qmGrpIDMap[(*it).first]; 01443 qmCustPCSizes[qmIndx] = (*it).second.size(); 01444 } 01445 01446 qmCustomPCIdxs = new int[qmTotCustPCs]; 01447 01448 if (simParams->qmCustomPCSel) { 01449 01450 int auxIter = 0; 01451 for (int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) { 01452 01453 Real currQMID = qmGrpID[grpIndx]; 01454 01455 iout << iINFO << "Loading " << qmPCVecMap[currQMID].size() 01456 << " custom point charges to QM Group " << grpIndx 01457 << " (ID: " << currQMID << ")\n" << endi; 01458 01459 for (int pcIndx=0; pcIndx<qmPCVecMap[currQMID].size(); pcIndx++) { 01460 qmCustomPCIdxs[auxIter] = qmPCVecMap[currQMID][pcIndx] ; 01461 auxIter++; 01462 } 01463 } 01464 } 01465 01466 // Conditional SMD option 01467 qmcSMD = simParams->qmCSMD; 01468 if (qmcSMD) { 01469 read_qm_csdm_file(qmGrpIDMap); 01470 } 01471 01474 01475 if (simParams->qmElecEmbed) { 01476 // Modifies Atom charges for Electrostatic Embedding. 01477 // QM atoms cannot have charges in the standard location, to keep 01478 // NAMD from calculating electrostatic interactions between QM and MM atoms. 01479 // We handle electrostatics ourselves in ComputeQM.C and in special 01480 // modifications for PME. 01481 for (int i=0; i<qmNumQMAtoms; i++) { 01482 qmAtmChrg[i] = atoms[qmAtmIndx[i]].charge; 01483 atoms[qmAtmIndx[i]].charge = 0; 01484 } 01485 } 01486 01487 01488 if ( simParams->extraBondsOn) { 01489 // Lifted from Molecule::build_extra_bonds 01490 01491 StringList *file = cfgList->find("extraBondsFile"); 01492 01493 char err_msg[512]; 01494 int a1,a2; float k, ref; 01495 01496 for ( ; file; file = file->next ) { // loop over files 01497 FILE *f = fopen(file->data,"r"); 01498 // if ( ! f ) { 01499 // sprintf(err_msg, "UNABLE TO OPEN EXTRA BONDS FILE %s", file->data); 01500 // NAMD_err(err_msg); 01501 // } else { 01502 // iout << iINFO << "READING EXTRA BONDS FILE " << file->data <<"\n"<<endi; 01503 // } 01504 01505 while ( 1 ) { 01506 char buffer[512]; 01507 int ret_code; 01508 do { 01509 ret_code = NAMD_read_line(f, buffer); 01510 } while ( (ret_code==0) && (NAMD_blank_string(buffer)) ); 01511 if (ret_code!=0) break; 01512 01513 char type[512]; 01514 sscanf(buffer,"%s",type); 01515 01516 int badline = 0; 01517 if ( ! strncasecmp(type,"bond",4) ) { 01518 if ( sscanf(buffer, "%s %d %d %f %f %s", 01519 type, &a1, &a2, &k, &ref, err_msg) != 5 ) badline = 1; 01520 01521 // If an extra bond is defined between QM atoms, we make 01522 // note so that it wont be deleted when we delete bonded 01523 // interactions between QM atoms. 01524 if( qmAtomGroup[a1] > 0 && qmAtomGroup[a2]) { 01525 Bond tmp; 01526 tmp.bond_type = 0; 01527 tmp.atom1 = a1; tmp.atom2 = a2; 01528 qmExtraBonds.add(tmp); 01529 } 01530 01531 01532 }else if ( ! strncasecmp(type,"#",1) ) { 01533 continue; // comment 01534 } 01535 01536 } 01537 fclose(f); 01538 } 01539 } 01540 01541 return; 01542 01543 #endif // MEM_OPT_VERSION 01544 }
| void Molecule::print_atoms | ( | Parameters * | params | ) |
Definition at line 5168 of file Molecule.C.
References DebugM, endi(), and Parameters::get_vdw_params().
Referenced by NamdState::loadStructure().
05169 { 05170 #ifdef MEM_OPT_VERSION 05171 DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi); 05172 #else 05173 register int i; 05174 Real sigma; 05175 Real epsilon; 05176 Real sigma14; 05177 Real epsilon14; 05178 05179 DebugM(2,"ATOM LIST\n" \ 05180 << "******************************************\n" \ 05181 << "NUM NAME TYPE RES MASS CHARGE CHARGE FEP-CHARGE" \ 05182 << "SIGMA EPSILON SIGMA14 EPSILON14\n" \ 05183 << endi); 05184 05185 for (i=0; i<numAtoms; i++) 05186 { 05187 params->get_vdw_params(&sigma, &epsilon, &sigma14, &epsilon14, 05188 atoms[i].vdw_type); 05189 05190 DebugM(2,i+1 << " " << atomNames[i].atomname \ 05191 << " " << atomNames[i].atomtype << " " \ 05192 << atomNames[i].resname << " " << atoms[i].mass \ 05193 << " " << atoms[i].charge << " " << sigma \ 05194 << " " << epsilon << " " << sigma14 \ 05195 << " " << epsilon14 << "\n" \ 05196 << endi); 05197 } 05198 #endif 05199 }
| void Molecule::print_bonds | ( | Parameters * | params | ) |
Definition at line 5211 of file Molecule.C.
References DebugM, endi(), and Parameters::get_bond_params().
Referenced by NamdState::loadStructure().
05212 { 05213 #ifdef MEM_OPT_VERSION 05214 DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi); 05215 #else 05216 register int i; 05217 Real k; 05218 Real x0; 05219 05220 DebugM(2,"BOND LIST\n" << "********************************\n" \ 05221 << "ATOM1 ATOM2 TYPE1 TYPE2 k x0" \ 05222 << endi); 05223 05224 for (i=0; i<numBonds; i++) 05225 { 05226 params->get_bond_params(&k, &x0, bonds[i].bond_type); 05227 05228 DebugM(2,bonds[i].atom1+1 << " " \ 05229 << bonds[i].atom2+1 << " " \ 05230 << atomNames[bonds[i].atom1].atomtype << " " \ 05231 << atomNames[bonds[i].atom2].atomtype << " " << k \ 05232 << " " << x0 << endi); 05233 } 05234 05235 #endif 05236 }
| void Molecule::print_exclusions | ( | ) |
Definition at line 5248 of file Molecule.C.
References DebugM, and endi().
Referenced by NamdState::loadStructure().
05249 { 05250 #ifdef MEM_OPT_VERSION 05251 DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi); 05252 #else 05253 register int i; 05254 05255 DebugM(2,"EXPLICIT EXCLUSION LIST\n" \ 05256 << "********************************\n" \ 05257 << "ATOM1 ATOM2 " \ 05258 << endi); 05259 05260 for (i=0; i<numExclusions; i++) 05261 { 05262 DebugM(2,exclusions[i].atom1+1 << " " \ 05263 << exclusions[i].atom2+1 << endi); 05264 } 05265 #endif 05266 }
| void Molecule::print_go_params | ( | ) |
Definition at line 548 of file GoMolecule.C.
References DebugM, go_array, j, MAX_GO_CHAINS, and NumGoChains.
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 }
| void Molecule::print_go_sigmas | ( | ) |
Definition at line 1134 of file GoMolecule.C.
References DebugM, goSigmaIndices, goSigmas, j, numAtoms, and numGoAtoms.
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 }
| void Molecule::put_stir_startTheta | ( | Real | theta, | |
| int | atomnum | |||
| ) | const [inline] |
Definition at line 1295 of file Molecule.h.
| void Molecule::read_go_file | ( | char * | fname | ) |
Definition at line 113 of file GoMolecule.C.
References go_val::cutoff, DebugM, endi(), go_val::epsilon, go_val::epsilonRep, go_val::exp_a, go_val::exp_b, go_val::exp_rep, FALSE, go_array, go_indices, iout, iWARN(), j, MAX_GO_CHAINS, MAX_RESTRICTIONS, NAMD_blank_string(), NAMD_die(), NAMD_find_first_word(), NAMD_read_line(), NumGoChains, go_val::restrictions, go_val::sigmaRep, and TRUE.
Referenced by build_go_params().
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 }
| void Molecule::read_parm | ( | Ambertoppar * | ) |
| void Molecule::receive_GoMolecule | ( | MIStream * | msg | ) |
Definition at line 1744 of file GoMolecule.C.
References atomChainTypes, go_val::cutoff, go_val::epsilon, go_val::epsilonRep, go_val::exp_a, go_val::exp_b, go_val::exp_rep, MIStream::get(), go_array, go_indices, goCoordinates, SimParameters::goForcesOn, goIndxLJA, goIndxLJB, SimParameters::goMethod, goNumLJPair, goResidIndices, goResids, goSigmaIndices, goSigmaPairA, goSigmaPairB, goSigmas, goWithinCutoff, j, MAX_GO_CHAINS, MAX_RESTRICTIONS, NAMD_die(), numAtoms, numGoAtoms, NumGoChains, pointerToGoBeg, pointerToGoEnd, go_val::restrictions, and go_val::sigmaRep.
01744 { 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 }
| void Molecule::receive_Molecule | ( | MIStream * | msg | ) |
Definition at line 5623 of file Molecule.C.
References SimParameters::alchOn, atomNamePool, atomSigPool, consForce, consForceIndexes, SimParameters::consForceOn, consTorqueIndexes, SimParameters::consTorqueOn, consTorqueParams, SimParameters::constraintsOn, DebugM, SimParameters::excludeFromPressure, SimParameters::fixedAtomsOn, gA, MIStream::get(), Parameters::get_num_vdw_params(), GridforceGrid::get_total_grids(), ObjectArena< Type >::getNewArray(), giSigma1, giSigma2, gMu1, gMu2, SimParameters::goGroPair, gRepulsive, gromacsPair_type, indxGaussA, indxGaussB, indxLJA, indxLJB, is_drude_psf, is_lonepairs_psf, isBFactorValid, isOccupancyValid, SimParameters::langevinOn, SimParameters::LCPOOn, SimParameters::lesOn, maxHydrogenGroupSize, maxMigrationGroupSize, SimParameters::mgridforceOn, SimParameters::movDragOn, numAcceptors, numCalcAngles, numCalcBonds, numCalcCrossterms, numCalcDihedrals, numCalcExclusions, numCalcFullExclusions, numCalcImpropers, numCalcLJPair, numConsForce, numConsTorque, numConstraints, numDonors, numExPressureAtoms, numFepFinal, numFepInitial, numFixedRigidBonds, numGaussPair, numGridforceGrids, numGridforces, numHydrogenGroups, numMigrationGroups, numMovDrag, numRotDrag, numStirredAtoms, pairC12, pairC6, SimParameters::pairInteractionOn, pointerToGaussBeg, pointerToGaussEnd, pointerToLJBeg, pointerToLJEnd, SimParameters::qmForcesOn, SimParameters::rotDragOn, SimParameters::soluteScalingOn, ss_index, ss_num_vdw_params, ss_vdw_type, SimParameters::stirOn, SimParameters::tCoupleOn, and GridforceGrid::unpack_grid().
Referenced by Node::resendMolecule().
05623 { 05624 // Get the atom information 05625 msg->get(numAtoms); 05626 05627 #ifdef MEM_OPT_VERSION 05628 //in the memory optimized version, only the atom signatures are recved 05629 //from the master Node. --Chao Mei 05630 05631 msg->get(massPoolSize); 05632 if(atomMassPool) delete [] atomMassPool; 05633 atomMassPool = new Real[massPoolSize]; 05634 msg->get(massPoolSize, atomMassPool); 05635 05636 msg->get(chargePoolSize); 05637 if(atomChargePool) delete [] atomChargePool; 05638 atomChargePool = new Real[chargePoolSize]; 05639 msg->get(chargePoolSize, atomChargePool); 05640 05641 //get atoms' signatures 05642 msg->get(atomSigPoolSize); 05643 if(atomSigPool) delete [] atomSigPool; 05644 atomSigPool = new AtomSignature[atomSigPoolSize]; 05645 for(int i=0; i<atomSigPoolSize; i++) 05646 atomSigPool[i].unpack(msg); 05647 05648 //get exclusions' signatures 05649 msg->get(exclSigPoolSize); 05650 if(exclSigPool) delete [] exclSigPool; 05651 exclSigPool = new ExclusionSignature[exclSigPoolSize]; 05652 for(int i=0; i<exclSigPoolSize; i++) 05653 exclSigPool[i].unpack(msg); 05654 05655 msg->get(numHydrogenGroups); 05656 msg->get(maxHydrogenGroupSize); 05657 msg->get(numMigrationGroups); 05658 msg->get(maxMigrationGroupSize); 05659 msg->get(isOccupancyValid); 05660 msg->get(isBFactorValid); 05661 05662 //get names for atoms 05663 msg->get(atomNamePoolSize); 05664 atomNamePool = new char *[atomNamePoolSize]; 05665 for(int i=0; i<atomNamePoolSize;i++) { 05666 int len; 05667 msg->get(len); 05668 atomNamePool[i] = nameArena->getNewArray(len+1); 05669 msg->get(len, atomNamePool[i]); 05670 } 05671 05672 if(simParams->fixedAtomsOn){ 05673 int numFixedAtomsSet; 05674 msg->get(numFixedAtoms); 05675 msg->get(numFixedAtomsSet); 05676 fixedAtomsSet = new AtomSetList(numFixedAtomsSet); 05677 msg->get(numFixedAtomsSet*sizeof(AtomSet), (char *)(fixedAtomsSet->begin())); 05678 } 05679 05680 if(simParams->constraintsOn){ 05681 int numConstrainedAtomsSet; 05682 msg->get(numConstraints); 05683 msg->get(numConstrainedAtomsSet); 05684 constrainedAtomsSet = new AtomSetList(numConstrainedAtomsSet); 05685 msg->get(numConstrainedAtomsSet*sizeof(AtomSet), (char *)(constrainedAtomsSet->begin())); 05686 } 05687 05688 #else 05689 delete [] atoms; 05690 atoms= new Atom[numAtoms]; 05691 msg->get(numAtoms*sizeof(Atom), (char*)atoms); 05692 05693 // Get the bond information 05694 msg->get(numRealBonds); 05695 msg->get(numBonds); 05696 if (numBonds) 05697 { 05698 delete [] bonds; 05699 bonds=new Bond[numBonds]; 05700 msg->get(numBonds*sizeof(Bond), (char*)bonds); 05701 } 05702 05703 // Get the angle information 05704 msg->get(numAngles); 05705 if (numAngles) 05706 { 05707 delete [] angles; 05708 angles=new Angle[numAngles]; 05709 msg->get(numAngles*sizeof(Angle), (char*)angles); 05710 } 05711 05712 // Get the dihedral information 05713 msg->get(numDihedrals); 05714 if (numDihedrals) 05715 { 05716 delete [] dihedrals; 05717 dihedrals=new Dihedral[numDihedrals]; 05718 msg->get(numDihedrals*sizeof(Dihedral), (char*)dihedrals); 05719 } 05720 05721 // Get the improper information 05722 msg->get(numImpropers); 05723 if (numImpropers) 05724 { 05725 delete [] impropers; 05726 impropers=new Improper[numImpropers]; 05727 msg->get(numImpropers*sizeof(Improper), (char*)impropers); 05728 } 05729 05730 // Get the crossterm information 05731 msg->get(numCrossterms); 05732 if (numCrossterms) 05733 { 05734 delete [] crossterms; 05735 crossterms=new Crossterm[numCrossterms]; 05736 msg->get(numCrossterms*sizeof(Crossterm), (char*)crossterms); 05737 } 05738 05739 // Get the hydrogen bond donors 05740 msg->get(numDonors); 05741 if (numDonors) 05742 { 05743 delete [] donors; 05744 donors=new Bond[numDonors]; 05745 msg->get(numDonors*sizeof(Bond), (char*)donors); 05746 } 05747 05748 // Get the hydrogen bond acceptors 05749 msg->get(numAcceptors); 05750 if (numAcceptors) 05751 { 05752 delete [] acceptors; 05753 acceptors=new Bond[numAcceptors]; 05754 msg->get(numAcceptors*sizeof(Bond), (char*)acceptors); 05755 } 05756 05757 // Get the exclusion information 05758 msg->get(numExclusions); 05759 if (numExclusions) 05760 { 05761 delete [] exclusions; 05762 exclusions=new Exclusion[numExclusions]; 05763 msg->get(numExclusions*sizeof(Exclusion), (char*)exclusions); 05764 } 05765 05766 // Get the constraint information, if they are active 05767 if (simParams->constraintsOn) 05768 { 05769 msg->get(numConstraints); 05770 05771 delete [] consIndexes; 05772 consIndexes = new int32[numAtoms]; 05773 05774 msg->get(numAtoms, consIndexes); 05775 05776 if (numConstraints) 05777 { 05778 delete [] consParams; 05779 consParams = new ConstraintParams[numConstraints]; 05780 05781 msg->get(numConstraints*sizeof(ConstraintParams), (char*)consParams); 05782 } 05783 } 05784 #endif 05785 05786 /* BEGIN gf */ 05787 if (simParams->mgridforceOn) 05788 { 05789 DebugM(3, "Receiving gridforce info\n"); 05790 05791 msg->get(numGridforceGrids); 05792 05793 DebugM(3, "numGridforceGrids = " << numGridforceGrids << "\n"); 05794 05795 delete [] numGridforces; 05796 numGridforces = new int[numGridforceGrids]; 05797 05798 delete [] gridfrcIndexes; // Should I be deleting elements of these first? 05799 delete [] gridfrcParams; 05800 delete [] gridfrcGrid; 05801 gridfrcIndexes = new int32*[numGridforceGrids]; 05802 gridfrcParams = new GridforceParams*[numGridforceGrids]; 05803 gridfrcGrid = new GridforceGrid*[numGridforceGrids]; 05804 05805 int grandTotalGrids = 0; 05806 for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) { 05807 msg->get(numGridforces[gridnum]); 05808 05809 gridfrcIndexes[gridnum] = new int32[numAtoms]; 05810 msg->get(numAtoms, gridfrcIndexes[gridnum]); 05811 05812 if (numGridforces[gridnum]) 05813 { 05814 gridfrcParams[gridnum] = new GridforceParams[numGridforces[gridnum]]; 05815 msg->get(numGridforces[gridnum]*sizeof(GridforceParams), (char*)gridfrcParams[gridnum]); 05816 } 05817 05818 gridfrcGrid[gridnum] = GridforceGrid::unpack_grid(gridnum, msg); 05819 05820 grandTotalGrids += gridfrcGrid[gridnum]->get_total_grids(); 05821 } 05822 } 05823 /* END gf */ 05824 05825 // Get the stirring information, if stirring is active 05826 if (simParams->stirOn) 05827 { 05828 msg->get(numStirredAtoms); 05829 05830 delete [] stirIndexes; 05831 stirIndexes = new int32[numAtoms]; 05832 05833 msg->get(numAtoms, stirIndexes); 05834 05835 if (numStirredAtoms) 05836 { 05837 delete [] stirParams; 05838 stirParams = new StirParams[numStirredAtoms]; 05839 05840 msg->get(numStirredAtoms*sizeof(StirParams), (char*)stirParams); 05841 } 05842 } 05843 05844 // Get the moving drag information, if it is active 05845 if (simParams->movDragOn) { 05846 msg->get(numMovDrag); 05847 delete [] movDragIndexes; 05848 movDragIndexes = new int32[numAtoms]; 05849 msg->get(numAtoms, movDragIndexes); 05850 if (numMovDrag) 05851 { 05852 delete [] movDragParams; 05853 movDragParams = new MovDragParams[numMovDrag]; 05854 msg->get(numMovDrag*sizeof(MovDragParams), (char*)movDragParams); 05855 } 05856 } 05857 05858 // Get the rotating drag information, if it is active 05859 if (simParams->rotDragOn) { 05860 msg->get(numRotDrag); 05861 delete [] rotDragIndexes; 05862 rotDragIndexes = new int32[numAtoms]; 05863 msg->get(numAtoms, rotDragIndexes); 05864 if (numRotDrag) 05865 { 05866 delete [] rotDragParams; 05867 rotDragParams = new RotDragParams[numRotDrag]; 05868 msg->get(numRotDrag*sizeof(RotDragParams), (char*)rotDragParams); 05869 } 05870 } 05871 05872 // Get the "constant" torque information, if it is active 05873 if (simParams->consTorqueOn) { 05874 msg->get(numConsTorque); 05875 delete [] consTorqueIndexes; 05876 consTorqueIndexes = new int32[numAtoms]; 05877 msg->get(numAtoms, consTorqueIndexes); 05878 if (numConsTorque) 05879 { 05880 delete [] consTorqueParams; 05881 consTorqueParams = new ConsTorqueParams[numConsTorque]; 05882 msg->get(numConsTorque*sizeof(ConsTorqueParams), (char*)consTorqueParams); 05883 } 05884 } 05885 05886 // Get the constant force information, if it's active 05887 if (simParams->consForceOn) 05888 { msg->get(numConsForce); 05889 delete [] consForceIndexes; 05890 consForceIndexes = new int32[numAtoms]; 05891 msg->get(numAtoms, consForceIndexes); 05892 if (numConsForce) 05893 { delete [] consForce; 05894 consForce = new Vector[numConsForce]; 05895 msg->get(numConsForce*sizeof(Vector), (char*)consForce); 05896 } 05897 } 05898 05899 if (simParams->excludeFromPressure) { 05900 exPressureAtomFlags = new int32[numAtoms]; 05901 msg->get(numExPressureAtoms); 05902 msg->get(numAtoms, exPressureAtomFlags); 05903 } 05904 05905 #ifndef MEM_OPT_VERSION 05906 // Get the langevin parameters, if they are active 05907 if (simParams->langevinOn || simParams->tCoupleOn) 05908 { 05909 delete [] langevinParams; 05910 langevinParams = new Real[numAtoms]; 05911 05912 msg->get(numAtoms, langevinParams); 05913 } 05914 05915 // Get the fixed atoms, if they are active 05916 if (simParams->fixedAtomsOn) 05917 { 05918 delete [] fixedAtomFlags; 05919 fixedAtomFlags = new int32[numAtoms]; 05920 05921 msg->get(numFixedAtoms); 05922 msg->get(numAtoms, fixedAtomFlags); 05923 msg->get(numFixedRigidBonds); 05924 } 05925 05926 if (simParams->qmForcesOn) 05927 { 05928 if( qmAtomGroup != 0) 05929 delete [] qmAtomGroup; 05930 qmAtomGroup = new Real[numAtoms]; 05931 05932 msg->get(numAtoms, qmAtomGroup); 05933 05934 msg->get(qmNumQMAtoms); 05935 05936 if( qmAtmChrg != 0) 05937 delete [] qmAtmChrg; 05938 qmAtmChrg = new Real[qmNumQMAtoms]; 05939 05940 msg->get(qmNumQMAtoms, qmAtmChrg); 05941 05942 if( qmAtmIndx != 0) 05943 delete [] qmAtmIndx; 05944 qmAtmIndx = new int[qmNumQMAtoms]; 05945 05946 msg->get(qmNumQMAtoms, qmAtmIndx); 05947 05948 msg->get(qmNoPC); 05949 05950 msg->get(qmNumBonds); 05951 05952 msg->get(qmMeNumBonds); 05953 05954 if( qmMeMMindx != 0) 05955 delete [] qmMeMMindx; 05956 qmMeMMindx = new int[qmMeNumBonds]; 05957 05958 msg->get(qmMeNumBonds, qmMeMMindx); 05959 05960 if( qmMeQMGrp != 0) 05961 delete [] qmMeQMGrp; 05962 qmMeQMGrp = new Real[qmMeNumBonds]; 05963 05964 msg->get(qmMeNumBonds, qmMeQMGrp); 05965 05966 msg->get(qmPCFreq); 05967 05968 msg->get(qmNumGrps); 05969 05970 if( qmGrpID != 0) 05971 delete [] qmGrpID; 05972 qmGrpID = new Real[qmNumGrps]; 05973 msg->get(qmNumGrps, qmGrpID); 05974 05975 if( qmCustPCSizes != 0) 05976 delete [] qmCustPCSizes; 05977 qmCustPCSizes = new int[qmNumGrps]; 05978 msg->get(qmNumGrps, qmCustPCSizes); 05979 05980 msg->get(qmTotCustPCs); 05981 05982 if( qmCustomPCIdxs != 0) 05983 delete [] qmCustomPCIdxs; 05984 qmCustomPCIdxs = new int[qmTotCustPCs]; 05985 msg->get(qmTotCustPCs, qmCustomPCIdxs); 05986 } 05987 05988 //fepb 05989 //receive fep atom info 05990 if (simParams->alchOn || simParams->lesOn || simParams->pairInteractionOn) { 05991 delete [] fepAtomFlags; 05992 fepAtomFlags = new unsigned char[numAtoms]; 05993 05994 msg->get(numFepInitial); 05995 msg->get(numFepFinal); 05996 msg->get(numAtoms*sizeof(unsigned char), (char*)fepAtomFlags); 05997 } 05998 //fepe 05999 06000 //soluteScaling 06001 if (simParams->soluteScalingOn) { 06002 delete [] ssAtomFlags; 06003 delete [] ss_vdw_type; 06004 delete [] ss_index; 06005 ssAtomFlags = new unsigned char[numAtoms]; 06006 ss_vdw_type = new int [params->get_num_vdw_params()]; 06007 ss_index = new int [numAtoms]; 06008 msg->get(numAtoms*sizeof(unsigned char), (char*)ssAtomFlags); 06009 msg->get(ss_num_vdw_params); 06010 msg->get(params->get_num_vdw_params()*sizeof(int), (char*)ss_vdw_type); 06011 msg->get(numAtoms*sizeof(int), (char*)ss_index); 06012 } 06013 //soluteScaling 06014 #ifdef OPENATOM_VERSION 06015 // This needs to be refactored into its own version 06016 if (simParams->openatomOn) { 06017 delete [] fepAtomFlags; 06018 fepAtomFlags = new unsigned char[numAtoms]; 06019 06020 msg->get(numFepInitial); 06021 msg->get(numAtoms*sizeof(unsigned char), (char*)fepAtomFlags); 06022 #endif //OPENATOM_VERSION 06023 06024 // DRUDE: receive data read from PSF 06025 msg->get(is_lonepairs_psf); 06026 if (is_lonepairs_psf) { 06027 msg->get(numLphosts); 06028 delete[] lphosts; 06029 lphosts = new Lphost[numLphosts]; 06030 msg->get(numLphosts*sizeof(Lphost), (char*)lphosts); 06031 } 06032 msg->get(is_drude_psf); 06033 if (is_drude_psf) { 06034 delete[] drudeConsts; 06035 drudeConsts = new DrudeConst[numAtoms]; 06036 msg->get(numAtoms*sizeof(DrudeConst), (char*)drudeConsts); 06037 msg->get(numAnisos); 06038 delete[] anisos; 06039 anisos = new Aniso[numAnisos]; 06040 msg->get(numAnisos*sizeof(Aniso), (char*)anisos); 06041 } 06042 // DRUDE 06043 06044 //LCPO 06045 if (simParams->LCPOOn) { 06046 delete [] lcpoParamType; 06047 lcpoParamType = new int[numAtoms]; 06048 msg->get(numAtoms, (int*)lcpoParamType); 06049 } 06050 06051 //Receive GromacsPairStuff -- JLai 06052 06053 if (simParams->goGroPair) { 06054 msg->get(numLJPair); 06055 delete [] indxLJA; 06056 indxLJA = new int[numLJPair]; 06057 msg->get(numLJPair,indxLJA); 06058 delete [] indxLJB; 06059 indxLJB = new int[numLJPair]; 06060 msg->get(numLJPair,indxLJB); 06061 delete [] pairC6; 06062 pairC6 = new Real[numLJPair]; 06063 msg->get(numLJPair,pairC6); 06064 delete [] pairC12; 06065 pairC12 = new Real[numLJPair]; 06066 msg->get(numLJPair,pairC12); 06067 delete [] gromacsPair_type; 06068 gromacsPair_type = new int[numLJPair]; 06069 msg->get(numLJPair,gromacsPair_type); 06070 delete [] pointerToLJBeg; 06071 pointerToLJBeg = new int[numAtoms]; 06072 msg->get((numAtoms),pointerToLJBeg); 06073 delete [] pointerToLJEnd; 06074 pointerToLJEnd = new int[numAtoms]; 06075 msg->get((numAtoms),pointerToLJEnd); 06076 // JLai 06077 delete [] gromacsPair; 06078 gromacsPair = new GromacsPair[numLJPair]; 06079 for(int i=0; i < numLJPair; i++) { 06080 gromacsPair[i].atom1 = indxLJA[i]; 06081 gromacsPair[i].atom2 = indxLJB[i]; 06082 gromacsPair[i].pairC6 = pairC6[i]; 06083 gromacsPair[i].pairC12 = pairC12[i]; 06084 gromacsPair[i].gromacsPair_type = gromacsPair_type[i]; 06085 } 06086 // 06087 msg->get(numGaussPair); 06088 delete [] indxGaussA; 06089 indxGaussA = new int[numGaussPair]; 06090 msg->get(numGaussPair,indxGaussA); 06091 delete [] indxGaussB; 06092 indxGaussB = new int[numGaussPair]; 06093 msg->get(numGaussPair,indxGaussB); 06094 delete [] gA; 06095 gA = new Real[numGaussPair]; 06096 msg->get(numGaussPair,gA); 06097 delete [] gMu1; 06098 gMu1 = new Real[numGaussPair]; 06099 msg->get(numGaussPair,gMu1); 06100 delete [] giSigma1; 06101 giSigma1 = new Real[numGaussPair]; 06102 msg->get(numGaussPair,giSigma1); 06103 delete [] gMu2; 06104 gMu2 = new Real[numGaussPair]; 06105 msg->get(numGaussPair,gMu2); 06106 delete [] giSigma2; 06107 giSigma2 = new Real[numGaussPair]; 06108 msg->get(numGaussPair,giSigma2); 06109 delete [] gRepulsive; 06110 gRepulsive = new Real[numGaussPair]; 06111 msg->get(numGaussPair,gRepulsive); 06112 delete [] pointerToGaussBeg; 06113 pointerToGaussBeg = new int[numAtoms]; 06114 msg->get((numAtoms),pointerToGaussBeg); 06115 delete [] pointerToGaussEnd; 06116 pointerToGaussEnd = new int[numAtoms]; 06117 msg->get((numAtoms),pointerToGaussEnd); 06118 // 06119 } 06120 #endif 06121 06122 // Now free the message 06123 delete msg; 06124 06125 #ifdef MEM_OPT_VERSION 06126 06127 build_excl_check_signatures(); 06128 06129 //set num{Calc}Tuples(Bonds,...,Impropers) to 0 06130 numBonds = numCalcBonds = 0; 06131 numAngles = numCalcAngles = 0; 06132 numDihedrals = numCalcDihedrals = 0; 06133 numImpropers = numCalcImpropers = 0; 06134 numCrossterms = numCalcCrossterms = 0; 06135 numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0; 06136 // JLai 06137 numLJPair = numCalcLJPair = 0; 06138 // End of JLai 06139 06140 #else 06141 06142 // analyze the data and find the status of each atom 06143 build_atom_status(); 06144 build_lists_by_atom(); 06145 06146 06147 #endif 06148 }
| void Molecule::reloadCharges | ( | float | charge[], | |
| int | n | |||
| ) |
Referenced by Node::reloadCharges().
| Real Molecule::rigid_bond_length | ( | int | atomnum | ) | const [inline] |
Definition at line 1440 of file Molecule.h.
Referenced by WorkDistrib::createAtomLists(), and outputCompressedFile().
| void Molecule::send_GoMolecule | ( | MOStream * | msg | ) |
Definition at line 1635 of file GoMolecule.C.
References atomChainTypes, go_val::cutoff, MOStream::end(), go_val::epsilon, go_val::epsilonRep, go_val::exp_a, go_val::exp_b, go_val::exp_rep, go_array, go_indices, goCoordinates, SimParameters::goForcesOn, goIndxLJA, goIndxLJB, SimParameters::goMethod, goNumLJPair, goResidIndices, goResids, goSigmaIndices, goSigmaPairA, goSigmaPairB, goSigmas, goWithinCutoff, j, MAX_GO_CHAINS, MAX_RESTRICTIONS, NAMD_die(), numAtoms, numGoAtoms, NumGoChains, pointerToGoBeg, pointerToGoEnd, MOStream::put(), go_val::restrictions, and go_val::sigmaRep.
01635 { 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 }
| void Molecule::send_Molecule | ( | MOStream * | msg | ) |
Definition at line 5279 of file Molecule.C.
References SimParameters::alchOn, atomNamePool, atomSigPool, consForce, consForceIndexes, SimParameters::consForceOn, consTorqueIndexes, SimParameters::consTorqueOn, consTorqueParams, SimParameters::constraintsOn, DebugM, MOStream::end(), endi(), SimParameters::excludeFromPressure, SimParameters::fixedAtomsOn, gA, Parameters::get_num_vdw_params(), giSigma1, giSigma2, gMu1, gMu2, SimParameters::goGroPair, gRepulsive, gromacsPair_type, indxGaussA, indxGaussB, indxLJA, indxLJB, is_drude_psf, is_lonepairs_psf, isBFactorValid, isOccupancyValid, SimParameters::langevinOn, SimParameters::LCPOOn, SimParameters::lesOn, maxHydrogenGroupSize, maxMigrationGroupSize, SimParameters::mgridforceOn, SimParameters::movDragOn, numAcceptors, numCalcAngles, numCalcBonds, numCalcCrossterms, numCalcDihedrals, numCalcExclusions, numCalcFullExclusions, numCalcImpropers, numCalcLJPair, numConsForce, numConsTorque, numConstraints, numDonors, numExPressureAtoms, numFepFinal, numFepInitial, numFixedRigidBonds, numGaussPair, numGridforceGrids, numGridforces, numHydrogenGroups, numMigrationGroups, numMovDrag, numRotDrag, numStirredAtoms, GridforceGrid::pack_grid(), pairC12, pairC6, SimParameters::pairInteractionOn, pointerToGaussBeg, pointerToGaussEnd, pointerToLJBeg, pointerToLJEnd, MOStream::put(), SimParameters::qmForcesOn, SimParameters::rotDragOn, SimParameters::soluteScalingOn, ss_index, ss_num_vdw_params, ss_vdw_type, SimParameters::stirOn, and SimParameters::tCoupleOn.
Referenced by Node::resendMolecule().
05279 { 05280 #ifdef MEM_OPT_VERSION 05281 //in the memory optimized version, only the atom signatures are broadcast 05282 //to other Nodes. --Chao Mei 05283 05284 msg->put(numAtoms); 05285 05286 msg->put(massPoolSize); 05287 msg->put(massPoolSize, atomMassPool); 05288 05289 msg->put(chargePoolSize); 05290 msg->put(chargePoolSize, atomChargePool); 05291 05292 //put atoms' signatures 05293 msg->put(atomSigPoolSize); 05294 for(int i=0; i<atomSigPoolSize; i++) 05295 atomSigPool[i].pack(msg); 05296 05297 //put atom's exclusion signatures 05298 msg->put(exclSigPoolSize); 05299 for(int i=0; i<exclSigPoolSize; i++) 05300 exclSigPool[i].pack(msg); 05301 05302 msg->put(numHydrogenGroups); 05303 msg->put(maxHydrogenGroupSize); 05304 msg->put(numMigrationGroups); 05305 msg->put(maxMigrationGroupSize); 05306 msg->put(isOccupancyValid); 05307 msg->put(isBFactorValid); 05308 05309 //put names for atoms 05310 msg->put(atomNamePoolSize); 05311 for(int i=0; i<atomNamePoolSize;i++) { 05312 int len = strlen(atomNamePool[i]); 05313 msg->put(len); 05314 msg->put(len*sizeof(char), atomNamePool[i]); 05315 } 05316 05317 if(simParams->fixedAtomsOn){ 05318 int numFixedAtomsSet = fixedAtomsSet->size(); 05319 msg->put(numFixedAtoms); 05320 msg->put(numFixedAtomsSet); 05321 msg->put(numFixedAtomsSet*sizeof(AtomSet), (char *)(fixedAtomsSet->begin())); 05322 } 05323 05324 if (simParams->constraintsOn) { 05325 int numConstrainedAtomsSet = constrainedAtomsSet->size(); 05326 msg->put(numConstraints); 05327 msg->put(numConstrainedAtomsSet); 05328 msg->put(numConstrainedAtomsSet*sizeof(AtomSet), (char *)(constrainedAtomsSet->begin())); 05329 } 05330 05331 #else 05332 msg->put(numAtoms); 05333 msg->put(numAtoms*sizeof(Atom), (char*)atoms); 05334 05335 // Send the bond information 05336 msg->put(numRealBonds); 05337 msg->put(numBonds); 05338 05339 if (numBonds) 05340 { 05341 msg->put(numBonds*sizeof(Bond), (char*)bonds); 05342 } 05343 05344 // Send the angle information 05345 msg->put(numAngles); 05346 if (numAngles) 05347 { 05348 msg->put(numAngles*sizeof(Angle), (char*)angles); 05349 } 05350 05351 // Send the dihedral information 05352 msg->put(numDihedrals); 05353 if (numDihedrals) 05354 { 05355 msg->put(numDihedrals*sizeof(Dihedral), (char*)dihedrals); 05356 } 05357 05358 // Send the improper information 05359 msg->put(numImpropers); 05360 if (numImpropers) 05361 { 05362 msg->put(numImpropers*sizeof(Improper), (char*)impropers); 05363 } 05364 05365 // Send the crossterm information 05366 msg->put(numCrossterms); 05367 if (numCrossterms) 05368 { 05369 msg->put(numCrossterms*sizeof(Crossterm), (char*)crossterms); 05370 } 05371 05372 // send the hydrogen bond donor information 05373 msg->put(numDonors); 05374 if(numDonors) 05375 { 05376 msg->put(numDonors*sizeof(Bond), (char*)donors); 05377 } 05378 05379 // send the hydrogen bond acceptor information 05380 msg->put(numAcceptors); 05381 if(numAcceptors) 05382 { 05383 msg->put(numAcceptors*sizeof(Bond), (char*)acceptors); 05384 } 05385 05386 // Send the exclusion information 05387 msg->put(numExclusions); 05388 if (numExclusions) 05389 { 05390 msg->put(numExclusions*sizeof(Exclusion), (char*)exclusions); 05391 } 05392 // Send the constraint information, if used 05393 if (simParams->constraintsOn) 05394 { 05395 msg->put(numConstraints); 05396 05397 msg->put(numAtoms, consIndexes); 05398 05399 if (numConstraints) 05400 { 05401 msg->put(numConstraints*sizeof(ConstraintParams), (char*)consParams); 05402 } 05403 } 05404 #endif 05405 05406 /* BEGIN gf */ 05407 // Send the gridforce information, if used 05408 if (simParams->mgridforceOn) 05409 { 05410 DebugM(3, "Sending gridforce info\n" << endi); 05411 msg->put(numGridforceGrids); 05412 05413 for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) { 05414 msg->put(numGridforces[gridnum]); 05415 msg->put(numAtoms, gridfrcIndexes[gridnum]); 05416 if (numGridforces[gridnum]) 05417 { 05418 msg->put(numGridforces[gridnum]*sizeof(GridforceParams), (char*)gridfrcParams[gridnum]); 05419 } 05420 GridforceGrid::pack_grid(gridfrcGrid[gridnum], msg); 05421 } 05422 } 05423 /* END gf */ 05424 05425 // Send the stirring information, if used 05426 if (simParams->stirOn) 05427 { 05428 //CkPrintf ("DEBUG: putting numStirredAtoms..\n"); 05429 msg->put(numStirredAtoms); 05430 //CkPrintf ("DEBUG: putting numAtoms,stirIndexes.. numAtoms=%d\n",numStirredAtoms); 05431 msg->put(numAtoms, stirIndexes); 05432 //CkPrintf ("DEBUG: if numStirredAtoms..\n"); 05433 if (numStirredAtoms) 05434 { 05435 //CkPrintf ("DEBUG: big put, with (char*)stirParams\n"); 05436 msg->put(numStirredAtoms*sizeof(StirParams), (char*)stirParams); 05437 } 05438 } 05439 05440 05441 // Send the moving drag information, if used 05442 if (simParams->movDragOn) { 05443 msg->put(numMovDrag); 05444 msg->put(numAtoms, movDragIndexes); 05445 if (numMovDrag) 05446 { 05447 msg->put(numMovDrag*sizeof(MovDragParams), (char*)movDragParams); 05448 } 05449 } 05450 05451 // Send the rotating drag information, if used 05452 if (simParams->rotDragOn) { 05453 msg->put(numRotDrag); 05454 msg->put(numAtoms, rotDragIndexes); 05455 if (numRotDrag) 05456 { 05457 msg->put(numRotDrag*sizeof(RotDragParams), (char*)rotDragParams); 05458 } 05459 } 05460 05461 // Send the "constant" torque information, if used 05462 if (simParams->consTorqueOn) { 05463 msg->put(numConsTorque); 05464 msg->put(numAtoms, consTorqueIndexes); 05465 if (numConsTorque) 05466 { 05467 msg->put(numConsTorque*sizeof(ConsTorqueParams), (char*)consTorqueParams); 05468 } 05469 } 05470 05471 // Send the constant force information, if used 05472 if (simParams->consForceOn) 05473 { msg->put(numConsForce); 05474 msg->put(numAtoms, consForceIndexes); 05475 if (numConsForce) 05476 msg->put(numConsForce*sizeof(Vector), (char*)consForce); 05477 } 05478 05479 if (simParams->excludeFromPressure) { 05480 msg->put(numExPressureAtoms); 05481 msg->put(numAtoms, exPressureAtomFlags); 05482 } 05483 05484 #ifndef MEM_OPT_VERSION 05485 // Send the langevin parameters, if active 05486 if (simParams->langevinOn || simParams->tCoupleOn) 05487 { 05488 msg->put(numAtoms, langevinParams); 05489 } 05490 05491 // Send fixed atoms, if active 05492 if (simParams->fixedAtomsOn) 05493 { 05494 msg->put(numFixedAtoms); 05495 msg->put(numAtoms, fixedAtomFlags); 05496 msg->put(numFixedRigidBonds); 05497 } 05498 05499 if (simParams->qmForcesOn) 05500 { 05501 msg->put(numAtoms, qmAtomGroup); 05502 msg->put(qmNumQMAtoms); 05503 msg->put(qmNumQMAtoms, qmAtmChrg); 05504 msg->put(qmNumQMAtoms, qmAtmIndx); 05505 msg->put(qmNoPC); 05506 msg->put(qmNumBonds); 05507 msg->put(qmMeNumBonds); 05508 msg->put(qmMeNumBonds, qmMeMMindx); 05509 msg->put(qmMeNumBonds, qmMeQMGrp); 05510 msg->put(qmPCFreq); 05511 msg->put(qmNumGrps); 05512 msg->put(qmNumGrps, qmGrpID); 05513 msg->put(qmNumGrps, qmCustPCSizes); 05514 msg->put(qmTotCustPCs); 05515 msg->put(qmTotCustPCs, qmCustomPCIdxs); 05516 } 05517 05518 //fepb 05519 // send fep atom info 05520 if (simParams->alchOn || simParams->lesOn || simParams->pairInteractionOn) { 05521 msg->put(numFepInitial); 05522 msg->put(numFepFinal); 05523 msg->put(numAtoms*sizeof(char), (char*)fepAtomFlags); 05524 } 05525 //fepe 05526 05527 if (simParams->soluteScalingOn) { 05528 msg->put(numAtoms*sizeof(char), (char*)ssAtomFlags); 05529 msg->put(ss_num_vdw_params); 05530 msg->put(params->get_num_vdw_params()*sizeof(int), (char*)ss_vdw_type); 05531 msg->put(numAtoms*sizeof(int), (char*)ss_index); 05532 } 05533 05534 #ifdef OPENATOM_VERSION 05535 // needs to be refactored into its own openatom version 05536 if (simParams->openatomOn ) { 05537 msg->put(numFepInitial); 05538 msg->put(numAtoms*sizeof(char), (char*)fepAtomFlags); 05539 } 05540 #endif //OPENATOM_VERSION 05541 05542 // DRUDE: send data read from PSF 05543 msg->put(is_lonepairs_psf); 05544 if (is_lonepairs_psf) { 05545 msg->put(numLphosts); 05546 msg->put(numLphosts*sizeof(Lphost), (char*)lphosts); 05547 } 05548 msg->put(is_drude_psf); 05549 if (is_drude_psf) { 05550 msg->put(numAtoms*sizeof(DrudeConst), (char*)drudeConsts); 05551 msg->put(numAnisos); 05552 msg->put(numAnisos*sizeof(Aniso), (char*)anisos); 05553 } 05554 // DRUDE 05555 05556 //LCPO 05557 if (simParams->LCPOOn) { 05558 msg->put(numAtoms, (int*)lcpoParamType); 05559 } 05560 05561 //Send GromacsPairStuff -- JLai 05562 if (simParams->goGroPair) { 05563 msg->put(numLJPair); 05564 msg->put(numLJPair,indxLJA); 05565 msg->put(numLJPair,indxLJB); 05566 msg->put(numLJPair,pairC6); 05567 msg->put(numLJPair,pairC12); 05568 msg->put(numLJPair,gromacsPair_type); 05569 msg->put((numAtoms),pointerToLJBeg); 05570 msg->put((numAtoms),pointerToLJEnd); 05571 msg->put(numGaussPair); 05572 msg->put(numGaussPair,indxGaussA); 05573 msg->put(numGaussPair,indxGaussB); 05574 msg->put(numGaussPair,gA); 05575 msg->put(numGaussPair,gMu1); 05576 msg->put(numGaussPair,giSigma1); 05577 msg->put(numGaussPair,gMu2); 05578 msg->put(numGaussPair,giSigma2); 05579 msg->put(numGaussPair,gRepulsive); 05580 msg->put((numAtoms),pointerToGaussBeg); 05581 msg->put((numAtoms),pointerToGaussEnd); 05582 } 05583 #endif 05584 05585 // Broadcast the message to the other nodes 05586 msg->end(); 05587 delete msg; 05588 05589 #ifdef MEM_OPT_VERSION 05590 05591 build_excl_check_signatures(); 05592 05593 //set num{Calc}Tuples(Bonds,...,Impropers) to 0 05594 numBonds = numCalcBonds = 0; 05595 numAngles = numCalcAngles = 0; 05596 numDihedrals = numCalcDihedrals = 0; 05597 numImpropers = numCalcImpropers = 0; 05598 numCrossterms = numCalcCrossterms = 0; 05599 numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0; 05600 // JLai 05601 numLJPair = numCalcLJPair = 0; 05602 // End of JLai 05603 05604 #else 05605 05606 // Now build arrays of indexes into these arrays by atom 05607 build_lists_by_atom(); 05608 05609 #endif 05610 }
| int Molecule::set_gridfrc_grid | ( | int | gridnum, | |
| GridforceGrid * | grid | |||
| ) | [inline] |
Definition at line 1272 of file Molecule.h.
References numGridforceGrids.
Referenced by Node::reloadGridforceGrid().
01273 { 01274 if (grid && gridnum >= 0 && gridnum < numGridforceGrids) { 01275 gridfrcGrid[gridnum] = grid; 01276 return 0; 01277 } else { 01278 return -1; 01279 } 01280 }
| void Molecule::set_qm_replaceAll | ( | Bool | newReplaceAll | ) | [inline] |
Definition at line 809 of file Molecule.h.
Referenced by NamdState::loadStructure().
| void Molecule::setBFactorData | ( | molfile_atom_t * | atomarray | ) |
Definition at line 3074 of file Molecule.C.
Referenced by Molecule().
| void Molecule::setOccupancyData | ( | molfile_atom_t * | atomarray | ) |
Definition at line 3067 of file Molecule.C.
Referenced by Molecule().
friend class AngleElem [friend] |
Definition at line 215 of file Molecule.h.
friend class AnisoElem [friend] |
Definition at line 219 of file Molecule.h.
friend class BondElem [friend] |
Definition at line 214 of file Molecule.h.
friend class CrosstermElem [friend] |
Definition at line 220 of file Molecule.h.
friend class DihedralElem [friend] |
Definition at line 216 of file Molecule.h.
friend class ExclElem [friend] |
Definition at line 213 of file Molecule.h.
friend class GromacsPairElem [friend] |
Definition at line 222 of file Molecule.h.
friend class ImproperElem [friend] |
Definition at line 217 of file Molecule.h.
friend class TholeElem [friend] |
Definition at line 218 of file Molecule.h.
friend class WorkDistrib [friend] |
Definition at line 224 of file Molecule.h.
Definition at line 561 of file Molecule.h.
Referenced by NamdState::loadStructure().
Definition at line 562 of file Molecule.h.
Referenced by NamdState::loadStructure().
Definition at line 563 of file Molecule.h.
Referenced by NamdState::loadStructure().
Definition at line 635 of file Molecule.h.
Referenced by build_go_arrays(), build_go_sigmas(), build_go_sigmas2(), get_go_force(), get_go_force2(), get_go_force_new(), goInit(), receive_GoMolecule(), and send_GoMolecule().
Definition at line 605 of file Molecule.h.
Referenced by ComputeConsForce::doForce(), receive_Molecule(), ComputeMgr::recvComputeConsForceMsg(), and send_Molecule().
Definition at line 604 of file Molecule.h.
Referenced by ComputeConsForce::doForce(), receive_Molecule(), ComputeMgr::recvComputeConsForceMsg(), and send_Molecule().
Definition at line 607 of file Molecule.h.
Referenced by ComputeConsTorque::doForce(), get_constorque_params(), is_atom_constorqued(), receive_Molecule(), and send_Molecule().
| ConsTorqueParams* Molecule::consTorqueParams |
Definition at line 608 of file Molecule.h.
Referenced by get_constorque_params(), receive_Molecule(), and send_Molecule().
Definition at line 676 of file Molecule.h.
Referenced by build_go_arrays(), and goInit().
Definition at line 677 of file Molecule.h.
Referenced by build_go_arrays(), and goInit().
Definition at line 668 of file Molecule.h.
Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().
Definition at line 670 of file Molecule.h.
Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().
Definition at line 672 of file Molecule.h.
Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().
Definition at line 669 of file Molecule.h.
Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().
Definition at line 671 of file Molecule.h.
Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().
| GoValue Molecule::go_array[MAX_GO_CHAINS *MAX_GO_CHAINS] |
Definition at line 1571 of file Molecule.h.
Referenced by get_go_cutoff(), get_go_epsilon(), get_go_epsilonRep(), get_go_exp_a(), get_go_exp_b(), get_go_exp_rep(), get_go_sigmaRep(), go_restricted(), print_go_params(), read_go_file(), receive_GoMolecule(), and send_GoMolecule().
| int Molecule::go_indices[MAX_GO_CHAINS+1] |
Definition at line 1572 of file Molecule.h.
Referenced by read_go_file(), receive_GoMolecule(), and send_GoMolecule().
Definition at line 640 of file Molecule.h.
Referenced by build_go_arrays(), get_go_force_new(), goInit(), receive_GoMolecule(), and send_GoMolecule().
| int* Molecule::goIndxLJA |
Definition at line 645 of file Molecule.h.
Referenced by build_go_sigmas2(), receive_GoMolecule(), and send_GoMolecule().
| int* Molecule::goIndxLJB |
Definition at line 646 of file Molecule.h.
Referenced by build_go_sigmas2(), get_go_force2(), receive_GoMolecule(), and send_GoMolecule().
Definition at line 644 of file Molecule.h.
Referenced by build_go_sigmas2(), receive_GoMolecule(), and send_GoMolecule().
Definition at line 642 of file Molecule.h.
Referenced by build_go_arrays(), build_go_sigmas(), build_go_sigmas2(), and goInit().
Definition at line 637 of file Molecule.h.
Referenced by build_go_sigmas2(), get_go_force2(), receive_GoMolecule(), and send_GoMolecule().
| int* Molecule::goResids |
Definition at line 641 of file Molecule.h.
Referenced by build_go_arrays(), get_go_force_new(), goInit(), receive_GoMolecule(), and send_GoMolecule().
Definition at line 636 of file Molecule.h.
Referenced by build_go_arrays(), build_go_sigmas(), build_go_sigmas2(), get_go_force(), get_go_force_new(), goInit(), print_go_sigmas(), receive_GoMolecule(), and send_GoMolecule().
| double* Molecule::goSigmaPairA |
Definition at line 647 of file Molecule.h.
Referenced by build_go_sigmas2(), get_go_force2(), receive_GoMolecule(), and send_GoMolecule().
| double* Molecule::goSigmaPairB |
Definition at line 648 of file Molecule.h.
Referenced by build_go_sigmas2(), get_go_force2(), receive_GoMolecule(), and send_GoMolecule().
Definition at line 638 of file Molecule.h.
Referenced by build_go_sigmas(), get_go_force(), goInit(), print_go_sigmas(), receive_GoMolecule(), and send_GoMolecule().
| bool* Molecule::goWithinCutoff |
Definition at line 639 of file Molecule.h.
Referenced by build_go_sigmas(), get_go_force(), goInit(), receive_GoMolecule(), and send_GoMolecule().
Definition at line 673 of file Molecule.h.
Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().
Definition at line 661 of file Molecule.h.
Referenced by receive_Molecule(), and send_Molecule().
Definition at line 631 of file Molecule.h.
Referenced by WorkDistrib::createAtomLists(), and outputCompressedFile().
| int* Molecule::indxGaussA |
Definition at line 666 of file Molecule.h.
Referenced by receive_Molecule(), and send_Molecule().
| int* Molecule::indxGaussB |
Definition at line 667 of file Molecule.h.
Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().
| int* Molecule::indxLJA |
Definition at line 657 of file Molecule.h.
Referenced by receive_Molecule(), and send_Molecule().
| int* Molecule::indxLJB |
Definition at line 658 of file Molecule.h.
Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().
Definition at line 459 of file Molecule.h.
Referenced by receive_Molecule(), and send_Molecule().
Definition at line 460 of file Molecule.h.
Referenced by receive_Molecule(), and send_Molecule().
Definition at line 1452 of file Molecule.h.
Referenced by receive_Molecule(), and send_Molecule().
Definition at line 1452 of file Molecule.h.
Referenced by receive_Molecule(), and send_Molecule().
Definition at line 592 of file Molecule.h.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), NamdState::loadStructure(), outputCompressedFile(), receive_Molecule(), and send_Molecule().
Definition at line 594 of file Molecule.h.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), NamdState::loadStructure(), outputCompressedFile(), receive_Molecule(), and send_Molecule().
Definition at line 567 of file Molecule.h.
Referenced by receive_Molecule(), and send_Molecule().
Definition at line 558 of file Molecule.h.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), ParallelIOMgr::bcastMolInfo(), buildAngleData(), delete_qm_bonded(), dumpbench(), AngleElem::getMoleculePointers(), NamdState::loadStructure(), Molecule(), and ParallelIOMgr::recvMolInfo().
Number of anisotropic terms.
Definition at line 575 of file Molecule.h.
Referenced by AnisoElem::getMoleculePointers().
Definition at line 554 of file Molecule.h.
Referenced by GlobalMasterFreeEnergy::addForce(), GlobalMasterEasy::addForce(), WorkDistrib::assignNodeToPatch(), ParallelIOMgr::bcastHydroBasedCounter(), build12Excls(), build13Excls(), build14Excls(), ComputeNonbondedCUDA::build_exclusions(), build_go_arrays(), build_go_sigmas(), build_go_sigmas2(), buildAtomData(), buildBondData(), buildExclusions(), colvarproxy_namd::calculate(), colvarproxy_namd::check_atom_id(), Controller::compareChecksums(), ComputeMgr::createComputes(), dumpbench(), GlobalMasterFreeEnergy::getMass(), GlobalMasterEasy::getMass(), GlobalMasterSymmetry::GlobalMasterSymmetry(), GlobalMasterTMD::GlobalMasterTMD(), colvarproxy_namd::init_atom_group(), integrateAllAtomSigs(), loadMolInfo(), NamdState::loadStructure(), Molecule(), num_deg_freedom(), outputCompressedFile(), WorkDistrib::patchMapInit(), prepare_qm(), print_go_sigmas(), receive_GoMolecule(), Controller::receivePressure(), ParallelIOMgr::recvAtomsCntPerPatch(), ComputeMgr::recvComputeConsForceMsg(), ComputeMsmSerialMgr::recvCoord(), ComputeGBISserMgr::recvCoord(), ComputeFmmSerialMgr::recvCoord(), ComputeExtMgr::recvCoord(), ComputeQMMgr::recvPartQM(), Node::reloadCharges(), GlobalMasterFreeEnergy::requestAtom(), GlobalMasterEasy::requestAtom(), Node::resendMolecule(), Node::resendMolecule2(), send_GoMolecule(), Node::startup(), Tcl_centerOfMass(), Tcl_centerOfNumber(), Tcl_loadCoords(), Tcl_radiusOfGyration(), wrap_coor_int(), and ~Molecule().
Definition at line 557 of file Molecule.h.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), ParallelIOMgr::bcastMolInfo(), buildBondData(), delete_qm_bonded(), dumpbench(), BondElem::getMoleculePointers(), NamdState::loadStructure(), Molecule(), prepare_qm(), and ParallelIOMgr::recvMolInfo().
Definition at line 614 of file Molecule.h.
Referenced by ParallelIOMgr::bcastMolInfo(), Controller::compareChecksums(), dumpbench(), receive_Molecule(), ParallelIOMgr::recvMolInfo(), and send_Molecule().
Definition at line 623 of file Molecule.h.
Referenced by Controller::compareChecksums().
Definition at line 613 of file Molecule.h.
Referenced by ParallelIOMgr::bcastMolInfo(), Controller::compareChecksums(), dumpbench(), receive_Molecule(), ParallelIOMgr::recvMolInfo(), and send_Molecule().
Definition at line 617 of file Molecule.h.
Referenced by ParallelIOMgr::bcastMolInfo(), Controller::compareChecksums(), receive_Molecule(), ParallelIOMgr::recvMolInfo(), and send_Molecule().
Definition at line 615 of file Molecule.h.
Referenced by ParallelIOMgr::bcastMolInfo(), Controller::compareChecksums(), dumpbench(), receive_Molecule(), ParallelIOMgr::recvMolInfo(), and send_Molecule().
Definition at line 618 of file Molecule.h.
Referenced by ParallelIOMgr::bcastMolInfo(), Controller::compareChecksums(), dumpbench(), receive_Molecule(), and send_Molecule().
Definition at line 619 of file Molecule.h.
Referenced by ParallelIOMgr::bcastMolInfo(), Controller::compareChecksums(), receive_Molecule(), and send_Molecule().
Definition at line 616 of file Molecule.h.
Referenced by ParallelIOMgr::bcastMolInfo(), Controller::compareChecksums(), dumpbench(), receive_Molecule(), ParallelIOMgr::recvMolInfo(), and send_Molecule().
Definition at line 654 of file Molecule.h.
Referenced by receive_Molecule(), and send_Molecule().
Definition at line 622 of file Molecule.h.
Referenced by Controller::compareChecksums().
Definition at line 603 of file Molecule.h.
Referenced by NamdState::loadStructure(), receive_Molecule(), and send_Molecule().
Definition at line 587 of file Molecule.h.
Referenced by is_atom_constorqued(), receive_Molecule(), and send_Molecule().
Definition at line 580 of file Molecule.h.
Referenced by is_atom_constrained(), NamdState::loadStructure(), num_deg_freedom(), num_group_deg_freedom(), receive_Molecule(), Controller::receivePressure(), and send_Molecule().
Definition at line 565 of file Molecule.h.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), ParallelIOMgr::bcastMolInfo(), buildCrosstermData(), delete_qm_bonded(), CrosstermElem::getMoleculePointers(), NamdState::loadStructure(), Molecule(), and ParallelIOMgr::recvMolInfo().
Definition at line 559 of file Molecule.h.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), ParallelIOMgr::bcastMolInfo(), buildDihedralData(), delete_qm_bonded(), dumpbench(), DihedralElem::getMoleculePointers(), NamdState::loadStructure(), Molecule(), and ParallelIOMgr::recvMolInfo().
Definition at line 566 of file Molecule.h.
Referenced by receive_Molecule(), and send_Molecule().
Number of Drude particles.
Definition at line 573 of file Molecule.h.
Referenced by NamdState::loadStructure(), and Controller::receivePressure().
Definition at line 568 of file Molecule.h.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), ExclElem::getMoleculePointers(), and NamdState::loadStructure().
Definition at line 590 of file Molecule.h.
Referenced by is_atom_exPressure(), receive_Molecule(), and send_Molecule().
Definition at line 600 of file Molecule.h.
Referenced by NamdState::loadStructure(), receive_Molecule(), and send_Molecule().
Definition at line 599 of file Molecule.h.
Referenced by NamdState::loadStructure(), num_deg_freedom(), receive_Molecule(), Controller::receivePressure(), and send_Molecule().
Definition at line 588 of file Molecule.h.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), ParallelIOMgr::integrateMigratedAtoms(), is_atom_fixed(), is_group_fixed(), NamdState::loadStructure(), num_fixed_atoms(), Controller::receivePressure(), and ParallelIOMgr::updateMolInfo().
Definition at line 595 of file Molecule.h.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), NamdState::loadStructure(), num_fixed_groups(), Controller::receivePressure(), and ParallelIOMgr::recvHydroBasedCounter().
Definition at line 597 of file Molecule.h.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), NamdState::loadStructure(), num_deg_freedom(), receive_Molecule(), Controller::receivePressure(), ParallelIOMgr::recvHydroBasedCounter(), and send_Molecule().
Definition at line 665 of file Molecule.h.
Referenced by receive_Molecule(), and send_Molecule().
Definition at line 634 of file Molecule.h.
Referenced by build_go_arrays(), build_go_sigmas(), build_go_sigmas2(), get_go_force(), goInit(), print_go_sigmas(), receive_GoMolecule(), and send_GoMolecule().
Definition at line 1573 of file Molecule.h.
Referenced by print_go_params(), read_go_file(), receive_GoMolecule(), and send_GoMolecule().
Definition at line 582 of file Molecule.h.
Referenced by build_gridforce_params(), ComputeGridForce::doForce(), get_gridfrc_grid(), is_atom_gridforced(), NamdState::loadStructure(), receive_Molecule(), send_Molecule(), and set_gridfrc_grid().
Definition at line 583 of file Molecule.h.
Referenced by build_gridforce_params(), receive_Molecule(), and send_Molecule().
Definition at line 591 of file Molecule.h.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), NamdState::loadStructure(), num_group_deg_freedom(), outputCompressedFile(), receive_Molecule(), Controller::receivePressure(), and send_Molecule().
Definition at line 564 of file Molecule.h.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), ParallelIOMgr::bcastMolInfo(), buildImproperData(), delete_qm_bonded(), dumpbench(), ImproperElem::getMoleculePointers(), NamdState::loadStructure(), Molecule(), and ParallelIOMgr::recvMolInfo().
Definition at line 653 of file Molecule.h.
Referenced by build_go_sigmas2(), and GromacsPairElem::getMoleculePointers().
Number of lone pairs.
Definition at line 572 of file Molecule.h.
Referenced by NamdState::loadStructure(), num_deg_freedom(), and Controller::receivePressure().
Number of lone pair host records in PSF.
Definition at line 576 of file Molecule.h.
Definition at line 593 of file Molecule.h.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), NamdState::loadStructure(), outputCompressedFile(), receive_Molecule(), and send_Molecule().
Definition at line 585 of file Molecule.h.
Referenced by is_atom_movdragged(), receive_Molecule(), and send_Molecule().
Definition at line 627 of file Molecule.h.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), and NamdState::loadStructure().
Definition at line 629 of file Molecule.h.
Referenced by ParallelIOMgr::bcastHydroBasedCounter(), and NamdState::loadStructure().
| int Molecule:: |