#include <Molecule.h>
|
||||||||||||
|
Definition at line 322 of file Molecule.C. References simParams. 00323 {
00324 initialize(simParams,param);
00325 }
|
|
||||||||||||||||||||
|
Definition at line 335 of file Molecule.C. References compress_psf_file(), SimParameters::genCompressedPsf, read_psf_file(), simParams, and SimParameters::useCompressedPsf. 00336 {
00337 initialize(simParams,param);
00338
00339 if(simParams->useCompressedPsf)
00340 read_compressed_psf_file(filename, param);
00341 else if(simParams->genCompressedPsf){
00342 compress_psf_file(this, filename, param, simParams, cfgList);
00343 }
00344 else
00345 read_psf_file(filename, param);
00346 }
|
|
||||||||||||||||
|
Definition at line 7550 of file Molecule.C. References Ambertoppar, read_parm(), and simParams. 07551 {
07552 initialize(simParams,param);
07553
07554 read_parm(amber_data);
07555 }
|
|
||||||||||||||||
|
Definition at line 7867 of file Molecule.C. References read_parm(), and simParams. 07869 {
07870 initialize(simParams,param);
07871
07872 read_parm(gromacsTopFile);
07873 }
|
|
|
Definition at line 360 of file Molecule.C. 00361 {
00362 /* Check to see if each array was ever allocated. If it was */
00363 /* then free it */
00364 if (atoms != NULL)
00365 delete [] atoms;
00366
00367 if (atomNames != NULL)
00368 {
00369 // subarrarys allocated from arena - automatically deleted
00370 delete [] atomNames;
00371 }
00372 delete nameArena;
00373
00374 if (resLookup != NULL)
00375 delete resLookup;
00376
00377 #ifdef MEM_OPT_VERSION
00378 if(eachAtomSig) delete [] eachAtomSig;
00379 if(atomSigPool) delete [] atomSigPool;
00380 #else
00381 if (bonds != NULL)
00382 delete [] bonds;
00383
00384 if (angles != NULL)
00385 delete [] angles;
00386
00387 if (dihedrals != NULL)
00388 delete [] dihedrals;
00389
00390 if (impropers != NULL)
00391 delete [] impropers;
00392
00393 if (crossterms != NULL)
00394 delete [] crossterms;
00395
00396 if (exclusions != NULL)
00397 delete [] exclusions;
00398 #endif
00399
00400 if (donors != NULL)
00401 delete [] donors;
00402
00403 if (acceptors != NULL)
00404 delete [] acceptors;
00405
00406 #ifdef MEM_OPT_VERSION
00407 if(exclSigPool) delete [] exclSigPool;
00408 if(exclChkSigPool) delete [] exclChkSigPool;
00409 if(eachAtomExclSig) delete [] eachAtomExclSig;
00410 #else
00411 if (bondsByAtom != NULL)
00412 delete [] bondsByAtom;
00413
00414 if (anglesByAtom != NULL)
00415 delete [] anglesByAtom;
00416
00417 if (dihedralsByAtom != NULL)
00418 delete [] dihedralsByAtom;
00419
00420 if (impropersByAtom != NULL)
00421 delete [] impropersByAtom;
00422
00423 if (crosstermsByAtom != NULL)
00424 delete [] crosstermsByAtom;
00425
00426 if (exclusionsByAtom != NULL)
00427 delete [] exclusionsByAtom;
00428
00429 if (fullExclusionsByAtom != NULL)
00430 delete [] fullExclusionsByAtom;
00431
00432 if (modExclusionsByAtom != NULL)
00433 delete [] modExclusionsByAtom;
00434
00435 if (all_exclusions != NULL)
00436 delete [] all_exclusions;
00437 #endif
00438
00439
00440 if (fixedAtomFlags != NULL)
00441 delete [] fixedAtomFlags;
00442
00443 if (stirIndexes != NULL)
00444 delete [] stirIndexes;
00445
00446
00447 #ifdef MEM_OPT_VERSION
00448 if(clusterSigs != NULL){
00449 delete [] clusterSigs;
00450 }
00451 #else
00452 if (cluster != NULL)
00453 delete [] cluster;
00454 if (clusterSize != NULL)
00455 delete [] clusterSize;
00456 #endif
00457
00458 if (exPressureAtomFlags != NULL)
00459 delete [] exPressureAtomFlags;
00460
00461 if (rigidBondLengths != NULL)
00462 delete [] rigidBondLengths;
00463
00464 //fepb
00465 if (fepAtomFlags != NULL)
00466 delete [] fepAtomFlags;
00467 //fepe
00468
00469
00470 #ifndef MEM_OPT_VERSION
00471 delete arena;
00472 delete exclArena;
00473 #endif
00474 }
|
|
|
Definition at line 461 of file Molecule.h. References atom_constants::charge, and Real. Referenced by NamdState::configListInit(), WorkDistrib::createAtomLists(), WorkDistrib::fillOnePatchAtoms(), print_atoms(), and Sequencer::reloadCharges(). 00462 {
00463 #ifdef MEM_OPT_VERSION
00464 return atomChargePool[eachAtomCharge[anum]];
00465 #else
00466 return(atoms[anum].charge);
00467 #endif
00468 }
|
|
|
Definition at line 451 of file Molecule.h. References atom_constants::mass, and Real. Referenced by NamdState::configListInit(), WorkDistrib::createAtomLists(), WorkDistrib::fillOnePatchAtoms(), GlobalMasterFreeEnergy::getMass(), GlobalMasterEasy::getMass(), print_atoms(), GlobalMaster::processData(), and ComputeGlobal::recvResults(). 00452 {
00453 #ifdef MEM_OPT_VERSION
00454 return atomMassPool[eachAtomMass[anum]];
00455 #else
00456 return(atoms[anum].mass);
00457 #endif
00458 }
|
|
|
Definition at line 471 of file Molecule.h. References Index, and atom_constants::vdw_type. Referenced by WorkDistrib::createAtomLists(), dumpbench(), WorkDistrib::fillOnePatchAtoms(), and SELF(). 00472 {
00473 return(atoms[anum].vdw_type);
00474 }
|
|
|
Definition at line 5855 of file Molecule.C. References PDB::atom(), consForce, consForceIndexes, int32, iout, iWARN(), NAMD_die(), PDB::num_atoms(), numConsForce, PDBAtom::occupancy(), Vector::x, PDBAtom::xcoor(), Vector::y, PDBAtom::ycoor(), Vector::z, and PDBAtom::zcoor(). Referenced by NamdState::configListInit(). 05856 { int i, index;
05857 PDB *forcePDB;
05858
05859 if (!filename) {
05860 // then all forces are zero to begin with; may be changed by
05861 // the consforceconfig command.
05862 iout << iWARN << "NO CONSTANT FORCES SPECIFIED, BUT CONSTANT FORCE IS ON . . .\n" << endi;
05863 consForceIndexes = new int32[numAtoms];
05864 for (i=0; i<numAtoms; i++) consForceIndexes[i] = -1;
05865 return;
05866 }
05867
05868 if ((forcePDB=new PDB(filename)) == NULL)
05869 NAMD_die("Memory allocation failed in Molecule::build_constant_forces");
05870 if (forcePDB->num_atoms() != numAtoms)
05871 NAMD_die("Number of atoms in constant force PDB doesn't match coordinate PDB");
05872
05873 // Allocate an array that will store an index into the constant force
05874 // array for each atom. If the atom has no constant force applied, its
05875 // value will be set to -1 in this array.
05876 consForceIndexes = new int32[numAtoms];
05877 if (consForceIndexes == NULL)
05878 NAMD_die("memory allocation failed in Molecule::build_constant_forces()");
05879
05880 // Loop through all the atoms and find out which ones have constant force
05881 numConsForce = 0;
05882 for (i=0; i<numAtoms; i++)
05883 if ((forcePDB->atom(i)->xcoor()==0 && forcePDB->atom(i)->ycoor()==0 &&
05884 forcePDB->atom(i)->zcoor()==0) || forcePDB->atom(i)->occupancy()==0)
05885 // This atom has no constant force
05886 consForceIndexes[i] = -1;
05887 else
05888 // This atom has constant force
05889 consForceIndexes[i] = numConsForce++;
05890
05891 if (numConsForce == 0)
05892 // Constant force was turned on, but there weren't really any non-zero forces
05893 iout << iWARN << "NO NON-ZERO FORCES WERE FOUND, BUT CONSTANT FORCE IS ON . . .\n" << endi;
05894 else
05895 { // Allocate an array to hold the forces
05896 consForce = new Vector[numConsForce];
05897 if (consForce == NULL)
05898 NAMD_die("memory allocation failed in Molecule::build_constant_forces");
05899 // Loop through all the atoms and assign the forces
05900 for (i=0; i<numAtoms; i++)
05901 if ((index=consForceIndexes[i]) != -1)
05902 { // This atom has constant force on it
05903 consForce[index].x = forcePDB->atom(i)->xcoor() * forcePDB->atom(i)->occupancy();
05904 consForce[index].y = forcePDB->atom(i)->ycoor() * forcePDB->atom(i)->occupancy();
05905 consForce[index].z = forcePDB->atom(i)->zcoor() * forcePDB->atom(i)->occupancy();
05906 }
05907 }
05908
05909 delete forcePDB;
05910 }
|
|
||||||||||||||||||||||||||||||||||||
|
Definition at line 5529 of file Molecule.C. References PDB::atom(), consTorqueIndexes, consTorqueParams, StringList::data, int32, iout, iWARN(), NAMD_die(), StringList::next, PDB::num_atoms(), numConsTorque, and Real. Referenced by NamdState::configListInit(). 05538 {
05539 PDB *tPDB, *aPDB, *pPDB, *vPDB; // Pointers to other PDB file(s)
05540 register int i; // Loop counter
05541 int current_index=0; // Index into values used
05542 int dtcol = 4; // Column to look for torque tag in
05543 Real dtval = 0; // Torque tag value retreived
05544 int dvcol = 4; // Column to look for angular velocity in
05545 Real dvval = 0; // Angular velocity value retreived
05546 char mainfilename[129]; // main "constant" torque PDB filename
05547 char axisfilename[129]; // "constant" torque axis PDB filename
05548 char pivotfilename[129]; // "constant" torque pivot point PDB filename
05549 char velfilename[129]; // "constant" torque angular velocity PDB filename
05550
05551 // Get the PDB to read the "constant" torque tags from. Again, if the
05552 // user gave us another file name, open that one. Otherwise, just
05553 // use the PDB with the initial coordinates
05554 if (consTorqueFile == NULL) {
05555 tPDB = initial_pdb;
05556
05557 } else {
05558
05559 if (consTorqueFile->next != NULL) {
05560 NAMD_die("Multiple definitions of \"constant\" torque tag file in configuration file");
05561 }
05562
05563 if ( (cwd == NULL) || (consTorqueFile->data[0] == '/') ) {
05564 strcpy(mainfilename, consTorqueFile->data);
05565 } else {
05566 strcpy(mainfilename, cwd);
05567 strcat(mainfilename, consTorqueFile->data);
05568 }
05569
05570 tPDB = new PDB(mainfilename);
05571 if ( tPDB == NULL ) {
05572 NAMD_die("Memory allocation failed in Molecule::build_constorque_params");
05573 }
05574
05575 if (tPDB->num_atoms() != numAtoms) {
05576 NAMD_die("Number of atoms in \"constant\" torque tag PDB doesn't match coordinate PDB");
05577 }
05578 }
05579
05580 // Get the PDB to read atom rotation axes. If no name given, use
05581 // consTorqueFile if both it AND consTorquePivotFile are defined. Can NOT
05582 // use the PDB coordinate file, nor consTorquePivotFile!
05583
05584 if (consTorqueAxisFile == NULL) {
05585 if (consTorqueFile == NULL) {
05586 NAMD_die("\"Constant\" torque axis file can not be same as coordinate PDB file");
05587 } else {
05588 if (consTorqueAxisFile->next != NULL) {
05589 NAMD_die("Multiple definitions of \"constant\" torque axis file in configuration file");
05590 };
05591 if (consTorquePivotFile == NULL) {
05592 NAMD_die("Need to specify at least one of consTorqueAxisFile and consTorquePivotFile; they can not be same");
05593 };
05594 aPDB = tPDB;
05595 };
05596
05597 } else {
05598
05599 if ( (cwd == NULL) || (consTorqueAxisFile->data[0] == '/') ) {
05600 strcpy(axisfilename, consTorqueAxisFile->data);
05601 } else {
05602 strcpy(axisfilename, cwd);
05603 strcat(axisfilename, consTorqueAxisFile->data);
05604 }
05605
05606 aPDB = new PDB(axisfilename);
05607 if ( aPDB == NULL ) {
05608 NAMD_die("Memory allocation failed in Molecule::build_constorque_params");
05609 }
05610
05611 if (aPDB->num_atoms() != numAtoms) {
05612 NAMD_die("Number of atoms in \"constant\" torque axis PDB doesn't match coordinate PDB");
05613 }
05614 };
05615
05616 // Get the PDB to read atom rotation pivot points. If no name given,
05617 // use consTorqueFile if both it AND consTorqueAxisFile are defined. Can
05618 // NOT use the PDB coordinate file, nor consTorqueAxisFile!
05619
05620 if (consTorquePivotFile == NULL) {
05621 if (consTorqueFile == NULL) {
05622 NAMD_die("\"Constant\" torque pivot point file can not be same as coordinate PDB file");
05623 } else {
05624 if (consTorquePivotFile->next != NULL) {
05625 NAMD_die("Multiple definitions of \"constant\" torque pivot point file in configuration file");
05626 };
05627 if (consTorqueAxisFile == NULL) {
05628 NAMD_die("Need to specify at least one of consTorqueAxisFile and consTorquePivotFile; they can not be same");
05629 };
05630 pPDB = tPDB;
05631 };
05632
05633 } else {
05634
05635 if ( (cwd == NULL) || (consTorquePivotFile->data[0] == '/') ) {
05636 strcpy(pivotfilename, consTorquePivotFile->data);
05637 } else {
05638 strcpy(pivotfilename, cwd);
05639 strcat(pivotfilename, consTorquePivotFile->data);
05640 }
05641
05642 pPDB = new PDB(pivotfilename);
05643 if ( pPDB == NULL ) {
05644 NAMD_die("Memory allocation failed in Molecule::build_constorque_params");
05645 }
05646
05647 if (pPDB->num_atoms() != numAtoms) {
05648 NAMD_die("Number of atoms in \"constant\" torque pivot point PDB doesn't match coordinate PDB");
05649 }
05650 };
05651
05652
05653 // Get the PDB to read atom angular velocities. If no name given,
05654 // use consTorqueFile (or the coordinate PDB file if consTorqueFile is not
05655 // defined).
05656
05657 if (consTorqueValFile == NULL) {
05658 vPDB = tPDB;
05659 } else {
05660 if (consTorqueValFile->next != NULL) {
05661 NAMD_die("Multiple definitions of \"constant\" torque velocity file in configuration file");
05662 };
05663
05664 if ( (cwd == NULL) || (consTorqueValFile->data[0] == '/') ) {
05665 strcpy(velfilename, consTorqueValFile->data);
05666 } else {
05667 strcpy(velfilename, cwd);
05668 strcat(velfilename, consTorqueValFile->data);
05669 }
05670
05671 vPDB = new PDB(velfilename);
05672 if ( vPDB == NULL ) {
05673 NAMD_die("Memory allocation failed in Molecule::build_constorque_params");
05674 }
05675
05676 if (vPDB->num_atoms() != numAtoms) {
05677 NAMD_die("Number of atoms in \"constant\" torque velocity PDB doesn't match coordinate PDB");
05678 }
05679 };
05680
05681 // Get the column that the torque tag is going to be in. If
05682 // consTorqueFile is defined, it can be in any of the 5 floating point
05683 // fields in the PDB (X, Y, Z, O, or B) which correspond to the
05684 // 1st, 2nd, ... 5th floating point fields. If consTorqueFile is NOT
05685 // defined, it can only be O or B fileds. The default is the O
05686 // (4th) field, which is the occupancy.
05687
05688 if (consTorqueCol == NULL) {
05689 dtcol = 4;
05690 } else {
05691 if (consTorqueCol->next != NULL) {
05692 NAMD_die("Multiple definitions of torque tag column in config file");
05693 };
05694
05695 if ( consTorqueFile == NULL
05696 && (!strcasecmp(consTorqueCol->data, "X")
05697 || !strcasecmp(consTorqueCol->data, "Y")
05698 || !strcasecmp(consTorqueCol->data, "Z"))) {
05699 NAMD_die("Can not read \"constant\" torque tags from X, Y, or Z column of the PDB coordinate file");
05700 };
05701 if (!strcasecmp(consTorqueCol->data, "X")) {
05702 dtcol=1;
05703 } else if (!strcasecmp(consTorqueCol->data, "Y")) {
05704 dtcol=2;
05705 } else if (!strcasecmp(consTorqueCol->data, "Z")) {
05706 dtcol=3;
05707 } else if (!strcasecmp(consTorqueCol->data, "O")) {
05708 dtcol=4;
05709 } else if (!strcasecmp(consTorqueCol->data, "B")) {
05710 dtcol=5;
05711 }
05712 else {
05713 NAMD_die("consTorqueCol must have value of X, Y, Z, O, or B");
05714 };
05715 };
05716
05717 // Get the column that the torque value is going to be
05718 // in. If consTorqueValFile is defined, it can be in any of the 5
05719 // floating point fields in the PDB (X, Y, Z, O, or B) which
05720 // correspond to the 1st, 2nd, ... 5th floating point fields. If
05721 // NEITHER of consTorqueValFile OR consTorqueFile is defined, it can
05722 // only be O or B fileds. The default is the O (4th) field, which
05723 // is the occupancy.
05724
05725 if (consTorqueValCol == NULL) {
05726 dvcol = 4;
05727 } else {
05728 if (consTorqueValCol->next != NULL) {
05729 NAMD_die("Multiple definitions of torque value column in config file");
05730 };
05731
05732 if (consTorqueValFile == NULL
05733 && consTorqueFile == NULL
05734 && strcasecmp(consTorqueCol->data, "B")
05735 && strcasecmp(consTorqueCol->data, "O")) {
05736 NAMD_die("Can not read \"constant\" torque values from X, Y, or Z column of the PDB coordinate file");
05737 };
05738 if (!strcasecmp(consTorqueValCol->data, "X")) {
05739 dvcol=1;
05740 } else if (!strcasecmp(consTorqueValCol->data, "Y")) {
05741 dvcol=2;
05742 } else if (!strcasecmp(consTorqueValCol->data, "Z")) {
05743 dvcol=3;
05744 } else if (!strcasecmp(consTorqueValCol->data, "O")) {
05745 dvcol=4;
05746 } else if (!strcasecmp(consTorqueValCol->data, "B")) {
05747 dvcol=5;
05748 }
05749 else {
05750 NAMD_die("consTorqueValCol must have value of X, Y, Z, O, or B");
05751 };
05752 };
05753
05754 // Allocate an array that will store an index into the torque
05755 // parameters for each atom. If the atom is not torqued, its
05756 // value will be set to -1 in this array.
05757 consTorqueIndexes = new int32[numAtoms];
05758 if (consTorqueIndexes == NULL) {
05759 NAMD_die("memory allocation failed in Molecule::build_constorque_params()");
05760 };
05761
05762 // Loop through all the atoms and find out which ones are torqued
05763 for (i=0; i<numAtoms; i++) {
05764 switch (dtcol) {
05765 case 1:
05766 dtval = (tPDB->atom(i))->xcoor();
05767 break;
05768 case 2:
05769 dtval = (tPDB->atom(i))->ycoor();
05770 break;
05771 case 3:
05772 dtval = (tPDB->atom(i))->zcoor();
05773 break;
05774 case 4:
05775 dtval = (tPDB->atom(i))->occupancy();
05776 break;
05777 case 5:
05778 dtval = (tPDB->atom(i))->temperaturefactor();
05779 break;
05780 }
05781
05782 if (dtval != 0.0) {
05783 // This atom is torqued
05784 consTorqueIndexes[i] = current_index;
05785 current_index++;
05786 } else {
05787 // This atom is not torqued
05788 consTorqueIndexes[i] = -1;
05789 }
05790 }
05791
05792 if (current_index == 0) {
05793 iout << iWARN << "NO TORQUED ATOMS WERE FOUND, BUT \"CONSTANT\" TORQUE IS ON . . . " << endi;
05794 } else {
05795 consTorqueParams = new ConsTorqueParams[current_index];
05796 if (consTorqueParams == NULL) {
05797 NAMD_die("memory allocation failed in Molecule::build_constorque_params");
05798 }
05799 };
05800
05801 numConsTorque = current_index;
05802
05803 // Loop through all the atoms and assign the parameters for those
05804 // that are torqued
05805 for (i=0; i<numAtoms; i++) {
05806 if (consTorqueIndexes[i] != -1) {
05807 consTorqueParams[consTorqueIndexes[i]].a[0] = (aPDB->atom(i))->xcoor();
05808 consTorqueParams[consTorqueIndexes[i]].a[1] = (aPDB->atom(i))->ycoor();
05809 consTorqueParams[consTorqueIndexes[i]].a[2] = (aPDB->atom(i))->zcoor();
05810 consTorqueParams[consTorqueIndexes[i]].p[0] = (pPDB->atom(i))->xcoor();
05811 consTorqueParams[consTorqueIndexes[i]].p[1] = (pPDB->atom(i))->ycoor();
05812 consTorqueParams[consTorqueIndexes[i]].p[2] = (pPDB->atom(i))->zcoor();
05813 switch (dvcol) {
05814 case 1:
05815 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->atom(i))->xcoor();
05816 break;
05817 case 2:
05818 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->atom(i))->ycoor();
05819 break;
05820 case 3:
05821 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->atom(i))->zcoor();
05822 break;
05823 case 4:
05824 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->atom(i))->occupancy();
05825 break;
05826 case 5:
05827 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->atom(i))->temperaturefactor();
05828 break;
05829 };
05830 };
05831 };
05832
05833 if (consTorqueFile != NULL) delete tPDB;
05834 if (consTorqueAxisFile != NULL) delete aPDB;
05835 if (consTorquePivotFile != NULL) delete pPDB;
05836 if (consTorqueValFile != NULL) delete vPDB;
05837 }
|
|
||||||||||||||||||||||||
|
Definition at line 4720 of file Molecule.C. References PDB::atom(), StringList::data, int32, iout, iWARN(), NAMD_die(), StringList::next, PDB::num_atoms(), numConstraints, and Real. Referenced by NamdState::configListInit(). 04726 {
04727 PDB *refPDB, *kPDB; // Pointer to other PDB's if used
04728 register int i; // Loop counter
04729 int current_index=0; // Index into values used
04730 int kcol = 4; // Column to look for force constant in
04731 Real kval = 0; // Force constant value retreived
04732 char filename[129]; // PDB filename
04733
04734 // Get the PDB object that contains the reference positions. If
04735 // the user gave another file name, use it. Otherwise, just use
04736 // the PDB file that has the initial coordinates. i.e., constrain
04737 // the atoms around their initial position. This is the most likely
04738 // case anyway
04739 if (consref == NULL)
04740 {
04741 refPDB = initial_pdb;
04742 }
04743 else
04744 {
04745 if (consref->next != NULL)
04746 {
04747 NAMD_die("Multiple definitions of constraint reference file in configruation file");
04748 }
04749
04750 if ( (cwd == NULL) || (consref->data[0] == '/') )
04751 {
04752 strcpy(filename, consref->data);
04753 }
04754 else
04755 {
04756 strcpy(filename, cwd);
04757 strcat(filename, consref->data);
04758 }
04759
04760 refPDB = new PDB(filename);
04761 if ( refPDB == NULL )
04762 {
04763 NAMD_die("Memory allocation failed in Molecule::build_constraint_params");
04764 }
04765
04766 if (refPDB->num_atoms() != numAtoms)
04767 {
04768 NAMD_die("Number of atoms in constraint reference PDB doesn't match coordinate PDB");
04769 }
04770 }
04771
04772 // Get the PDB to read the force constants from. Again, if the user
04773 // gave us another file name, open that one. Otherwise, just use
04774 // the PDB with the initial coordinates
04775 if (conskfile == NULL)
04776 {
04777 kPDB = initial_pdb;
04778 }
04779 else
04780 {
04781 if (conskfile->next != NULL)
04782 {
04783 NAMD_die("Multiple definitions of constraint constant file in configuration file");
04784 }
04785
04786 if ( (consref != NULL) && (strcasecmp(consref->data, conskfile->data) == 0) )
04787 {
04788 // Same PDB used for reference positions and force constants
04789 kPDB = refPDB;
04790 }
04791 else
04792 {
04793 if ( (cwd == NULL) || (conskfile->data[0] == '/') )
04794 {
04795 strcpy(filename, conskfile->data);
04796 }
04797 else
04798 {
04799 strcpy(filename, cwd);
04800 strcat(filename, conskfile->data);
04801 }
04802
04803 kPDB = new PDB(filename);
04804 if ( kPDB == NULL )
04805 {
04806 NAMD_die("Memory allocation failed in Molecule::build_constraint_params");
04807 }
04808
04809 if (kPDB->num_atoms() != numAtoms)
04810 {
04811 NAMD_die("Number of atoms in constraint constant PDB doesn't match coordinate PDB");
04812 }
04813 }
04814 }
04815
04816 // Get the column that the force constant is going to be in. It
04817 // can be in any of the 5 floating point fields in the PDB, according
04818 // to what the user wants. The allowable fields are X, Y, Z, O, or
04819 // B which correspond to the 1st, 2nd, ... 5th floating point fields.
04820 // The default is the 4th field, which is the occupancy
04821 if (conskcol == NULL)
04822 {
04823 kcol = 4;
04824 }
04825 else
04826 {
04827 if (conskcol->next != NULL)
04828 {
04829 NAMD_die("Multiple definitions of harmonic constraint column in config file");
04830 }
04831
04832 if (strcasecmp(conskcol->data, "X") == 0)
04833 {
04834 kcol=1;
04835 }
04836 else if (strcasecmp(conskcol->data, "Y") == 0)
04837 {
04838 kcol=2;
04839 }
04840 else if (strcasecmp(conskcol->data, "Z") == 0)
04841 {
04842 kcol=3;
04843 }
04844 else if (strcasecmp(conskcol->data, "O") == 0)
04845 {
04846 kcol=4;
04847 }
04848 else if (strcasecmp(conskcol->data, "B") == 0)
04849 {
04850 kcol=5;
04851 }
04852 else
04853 {
04854 NAMD_die("conskcol must have value of X, Y, Z, O, or B");
04855 }
04856 }
04857
04858 // Allocate an array that will store an index into the constraint
04859 // parameters for each atom. If the atom is not constrained, its
04860 // value will be set to -1 in this array.
04861 consIndexes = new int32[numAtoms];
04862
04863 if (consIndexes == NULL)
04864 {
04865 NAMD_die("memory allocation failed in Molecule::build_constraint_params()");
04866 }
04867
04868 // Loop through all the atoms and find out which ones are constrained
04869 for (i=0; i<numAtoms; i++)
04870 {
04871 // Get the k value based on where we were told to find it
04872 switch (kcol)
04873 {
04874 case 1:
04875 kval = (kPDB->atom(i))->xcoor();
04876 break;
04877 case 2:
04878 kval = (kPDB->atom(i))->ycoor();
04879 break;
04880 case 3:
04881 kval = (kPDB->atom(i))->zcoor();
04882 break;
04883 case 4:
04884 kval = (kPDB->atom(i))->occupancy();
04885 break;
04886 case 5:
04887 kval = (kPDB->atom(i))->temperaturefactor();
04888 break;
04889 }
04890
04891 if (kval > 0.0)
04892 {
04893 // This atom is constrained
04894 consIndexes[i] = current_index;
04895 current_index++;
04896 }
04897 else
04898 {
04899 // This atom is not constrained
04900 consIndexes[i] = -1;
04901 }
04902 }
04903
04904 if (current_index == 0)
04905 {
04906 // Constraints were turned on, but there weren't really any constrained
04907 iout << iWARN << "NO CONSTRAINED ATOMS WERE FOUND, BUT CONSTRAINTS ARE ON . . .\n" << endi;
04908 }
04909 else
04910 {
04911 // Allocate an array to hold the constraint parameters
04912 consParams = new ConstraintParams[current_index];
04913
04914 if (consParams == NULL)
04915 {
04916 NAMD_die("memory allocation failed in Molecule::build_constraint_params");
04917 }
04918 }
04919
04920 numConstraints = current_index;
04921
04922 // Loop through all the atoms and assign the parameters for those
04923 // that are constrained
04924 for (i=0; i<numAtoms; i++)
04925 {
04926 if (consIndexes[i] != -1)
04927 {
04928 // This atom is constrained, so get the k value again
04929 switch (kcol)
04930 {
04931 case 1:
04932 consParams[consIndexes[i]].k = (kPDB->atom(i))->xcoor();
04933 break;
04934 case 2:
04935 consParams[consIndexes[i]].k = (kPDB->atom(i))->ycoor();
04936 break;
04937 case 3:
04938 consParams[consIndexes[i]].k = (kPDB->atom(i))->zcoor();
04939 break;
04940 case 4:
04941 consParams[consIndexes[i]].k = (kPDB->atom(i))->occupancy();
04942 break;
04943 case 5:
04944 consParams[consIndexes[i]].k = (kPDB->atom(i))->temperaturefactor();
04945 break;
04946 }
04947
04948 // Get the reference position
04949 consParams[consIndexes[i]].refPos.x = (refPDB->atom(i))->xcoor();
04950 consParams[consIndexes[i]].refPos.y = (refPDB->atom(i))->ycoor();
04951 consParams[consIndexes[i]].refPos.z = (refPDB->atom(i))->zcoor();
04952 }
04953 }
04954
04955 // If we had to create new PDB objects, delete them now
04956 if (consref != NULL)
04957 {
04958 delete refPDB;
04959 }
04960
04961 if ((conskfile != NULL) &&
04962 !((consref != NULL) &&
04963 (strcasecmp(consref->data, conskfile->data) == 0)
04964 )
04965 )
04966 {
04967 delete kPDB;
04968 }
04969
04970 }
|
|
||||||||||||||||||||
|
Definition at line 6880 of file Molecule.C. References PDB::atom(), StringList::data, iINFO(), int32, iout, NAMD_die(), StringList::next, PDB::num_atoms(), numExPressureAtoms, and Real. Referenced by NamdState::configListInit(). 06881 {
06882
06883 PDB *bPDB; // Pointer to PDB object to use
06884 int bcol = 4; // Column that data is in
06885 Real bval = 0; // b value from PDB file
06886 int i; // Loop counter
06887 char filename[129]; // Filename
06888
06889 // Get the PDB object that contains the b values. If
06890 // the user gave another file name, use it. Otherwise, just use
06891 // the PDB file that has the initial coordinates.
06892 if (fixedfile == NULL) {
06893 bPDB = initial_pdb;
06894 } else {
06895 if (fixedfile->next != NULL) {
06896 NAMD_die("Multiple definitions of excluded pressure atoms PDB file in configuration file");
06897 }
06898
06899 if ( (cwd == NULL) || (fixedfile->data[0] == '/') ) {
06900 strcpy(filename, fixedfile->data);
06901 } else {
06902 strcpy(filename, cwd);
06903 strcat(filename, fixedfile->data);
06904 }
06905 bPDB = new PDB(filename);
06906 if ( bPDB == NULL ) {
06907 NAMD_die("Memory allocation failed in Molecule::build_exPressure_atoms");
06908 }
06909
06910 if (bPDB->num_atoms() != numAtoms) {
06911 NAMD_die("Number of atoms in excludedPressure atoms PDB doesn't match coordinate PDB");
06912 }
06913 }
06914
06915 // Get the column that the b vaules are in. It
06916 // can be in any of the 5 floating point fields in the PDB, according
06917 // to what the user wants. The allowable fields are X, Y, Z, O, or
06918 // B which correspond to the 1st, 2nd, ... 5th floating point fields.
06919 // The default is the 4th field, which is the occupancy
06920 if (fixedcol == NULL) {
06921 bcol = 4;
06922 } else {
06923 if (fixedcol->next != NULL) {
06924 NAMD_die("Multiple definitions of excludedPressure atoms column in config file");
06925 }
06926
06927 if (strcasecmp(fixedcol->data, "X") == 0) {
06928 bcol=1;
06929 } else if (strcasecmp(fixedcol->data, "Y") == 0) {
06930 bcol=2;
06931 } else if (strcasecmp(fixedcol->data, "Z") == 0) {
06932 bcol=3;
06933 } else if (strcasecmp(fixedcol->data, "O") == 0) {
06934 bcol=4;
06935 } else if (strcasecmp(fixedcol->data, "B") == 0) {
06936 bcol=5;
06937 } else {
06938 NAMD_die("excludedPressureFileCol must have value of X, Y, Z, O, or B");
06939 }
06940 }
06941
06942 // Allocate the array to hold all the data
06943 exPressureAtomFlags = new int32[numAtoms];
06944
06945 if (exPressureAtomFlags == NULL) {
06946 NAMD_die("memory allocation failed in Molecule::build_fixed_atoms()");
06947 }
06948
06949 numExPressureAtoms = 0;
06950
06951 // Loop through all the atoms and get the b value
06952 for (i=0; i<numAtoms; i++) {
06953 // Get the k value based on where we were told to find it
06954 switch (bcol) {
06955 case 1: bval = (bPDB->atom(i))->xcoor(); break;
06956 case 2: bval = (bPDB->atom(i))->ycoor(); break;
06957 case 3: bval = (bPDB->atom(i))->zcoor(); break;
06958 case 4: bval = (bPDB->atom(i))->occupancy(); break;
06959 case 5: bval = (bPDB->atom(i))->temperaturefactor(); break;
06960 }
06961
06962 // Assign the b value
06963 if ( bval != 0 ) {
06964 exPressureAtomFlags[i] = 1;
06965 numExPressureAtoms++;
06966 } else {
06967 exPressureAtomFlags[i] = 0;
06968 }
06969 }
06970 if (fixedfile != NULL)
06971 delete bPDB;
06972
06973 iout << iINFO << "Got " << numExPressureAtoms << " excluded pressure atoms."
06974 << endi;
06975 }
|
|
||||||||||||