--- a/tests/adhoc/TestMemoryUse.cpp +++ /dev/null @@ -1,170 +0,0 @@ -#include "Molmodel.h" -#include - -#if defined (__linux__) -# include - -#elif defined(__APPLE__) -# include -# include - -#elif defined(_WINDOWS) -# include -# include "psapi.h" - -#else -# include -#endif - -/// The amount of memory currently being used by this process, in bytes. -/// By default, returns the full virtual arena, but if resident=true, -/// it will report just the resident set in RAM (if supported on that OS). -size_t memory_used (bool resident=false) -{ -#if defined(__linux__) - // Ugh, getrusage doesn't work well on Linux. Try grabbing info - // directly from the /proc pseudo-filesystem. Reading from - // /proc/self/statm gives info on your own process, as one line of - // numbers that are: virtual mem program size, resident set size, - // shared pages, text/code, data/stack, library, dirty pages. The - // mem sizes should all be multiplied by the page size. - size_t size = 0; - FILE *file = fopen("/proc/self/statm", "r"); - if (file) { - unsigned long vm = 0; - fscanf (file, "%ul", &vm); // Just need the first num: vm size - fclose (file); - size = (size_t)vm * getpagesize(); - } - return size; - -#elif defined(__APPLE__) - // Inspired by: - // http://miknight.blogspot.com/2005/11/resident-set-size-in-mac-os-x.html - struct task_basic_info t_info; - mach_msg_type_number_t t_info_count = TASK_BASIC_INFO_COUNT; - task_info(current_task(), TASK_BASIC_INFO, (task_info_t)&t_info, &t_info_count); - size_t size = (resident ? t_info.resident_size : t_info.virtual_size); - return size; - -#elif defined(_WINDOWS) - // According to MSDN... - PROCESS_MEMORY_COUNTERS count; - if (GetProcessMemoryInfo (GetCurrentProcess(), &count, sizeof (count))) - return count.PagefileUsage; - else return 0; - -#else - // No idea what platform this is - return 0; // Punt -#endif -} - -#include -#include - -using namespace SimTK; -using namespace std; - -void testRnaResources(int sequenceLength) -{ - time_t initialTime; time(&initialTime); - size_t initialMemFootprint = memory_used(); - - char resBuf[100]; - RNA rna("", false); - for (int b = 0; b < sequenceLength; ++b) { - // itoa(b+1, resBuf, 10); - sprintf(resBuf, "%d", b+1); - rna.appendResidue(resBuf, RibonucleotideResidue::Adenylate().withPhosphodiester()); - } - - time_t compoundTime; time(&compoundTime); - size_t compoundMemFootprint = memory_used(); - - CompoundSystem system; - SimbodyMatterSubsystem matter(system); - GeneralForceSubsystem forces(system); - DuMMForceFieldSubsystem dumm(system); - dumm.loadAmber99Parameters(); - dumm.setAllGlobalScaleFactors(0); - dumm.setBondTorsionGlobalScaleFactor(1); - rna.assignBiotypes(); - system.adoptCompound(rna); - - time_t preModelTime; time(&preModelTime); - system.modelCompounds(); - time_t postModelTime; time(&postModelTime); - - State state = system.realizeTopology(); - system.realize(state, Stage::Position); - VerletIntegrator integrator(system, 1e-3); - TimeStepper timeStepper(system, integrator); - timeStepper.initialize(state); - timeStepper.stepTo(0.020); - - time_t simTime; time(&simTime); - size_t systemMemFootprint = memory_used(); - - printf("Compound: %.2f Mb/%d bases", (compoundMemFootprint - initialMemFootprint)/1e6, sequenceLength); - printf("\tSystem: %.2f Mb/%d bases\n", (systemMemFootprint - compoundMemFootprint)/1e6, sequenceLength); - - printf("Compound: %.2f s/%d bases", difftime(preModelTime, initialTime), sequenceLength); - printf("\tModel: %.2f s/%d bases", difftime(postModelTime, preModelTime), sequenceLength); - printf("\tSimulation: %.2f s/%d bases /0.020 ps\n", difftime(simTime, postModelTime), sequenceLength); -} - -void testMemoryUse() -{ - // testRnaResources(1); - testRnaResources(10); - // testRnaResources(100); - // testRnaResources(200); - testRnaResources(400); - // testRnaResources(800); - // testRnaResources(1600); - // testRnaResources(2000); - // testRnaResources(2200); - // testRnaResources(2300); - // testRnaResources(3000); - // testRnaResources(6000); - // testRnaResources(6200); - // testRnaResources(6250); - // testRnaResources(6260); - // testRnaResources(6270); - // testRnaResources(6280); - // testRnaResources(6282); - // testRnaResources(6284); - // testRnaResources(6285); - // testRnaResources(6286); - // testRnaResources(6288); - // testRnaResources(6290); - // testRnaResources(6300); - // testRnaResources(6350); - // testRnaResources(6400); - // testRnaResources(6450); - // testRnaResources(6500); - // testRnaResources(6700); - // testRnaResources(7000); - // testRnaResources(8000); - // testRnaResources(9000); - // testRnaResources(10000); - // testRnaResources(50000); -} - -int main() { - try { - testMemoryUse(); - cout << "PASSED" << endl; - return 0; - } - catch (const std::exception& e) - { - printf("EXCEPTION THROWN: %s\n", e.what()); - return 1; - } - catch (...) - { - printf("UNKNOWN EXCEPTION THROWN\n"); - } return 1; -} --- a/tests/TestDihedrals2.cpp +++ /dev/null @@ -1,111 +0,0 @@ - -#include "SimTKmolmodel.h" -#include "molmodel/internal/Pdb.h" - -#include - -using namespace SimTK; -using namespace std; - -int main() -{ - CompoundSystem system; - DuMMForceFieldSubsystem dumm(system); - dumm.loadAmber99Parameters(); - - const char* uudPdbString = "" -"ATOM 1 P G 1 -10.655 -7.494 9.276 1.00 0.00 P \n" -"ATOM 2 OP1 G 1 -11.450 -7.642 8.036 1.00 0.00 O \n" -"ATOM 3 OP2 G 1 -9.248 -7.952 9.308 1.00 0.00 O \n" -"ATOM 4 O5* G 1 -10.689 -5.944 9.713 1.00 0.00 O \n" -"ATOM 5 C5* G 1 -11.850 -5.146 9.440 1.00 0.00 C \n" -"ATOM 6 H5*1 G 1 -12.691 -5.536 9.997 1.00 0.00 H \n" -"ATOM 7 H5*2 G 1 -12.064 -5.200 8.363 1.00 0.00 H \n" -"ATOM 8 C4* G 1 -11.678 -3.694 9.863 1.00 0.00 C \n" -"ATOM 9 H4* G 1 -12.581 -3.135 9.623 1.00 0.00 H \n" -"ATOM 10 O4* G 1 -11.419 -3.590 11.266 1.00 0.00 O \n" -"ATOM 11 C3* G 1 -10.481 -3.050 9.181 1.00 0.00 C \n" -"ATOM 12 H3* G 1 -9.675 -3.764 9.008 1.00 0.00 H \n" -"ATOM 13 O3* G 1 -10.924 -2.423 7.974 1.00 0.00 O \n" -"ATOM 14 C2* G 1 -10.085 -1.969 10.167 1.00 0.00 C \n" -"ATOM 15 H2* G 1 -9.021 -1.742 10.071 1.00 0.00 H \n" -"ATOM 16 C1* G 1 -10.357 -2.648 11.503 1.00 0.00 C \n" -"ATOM 17 H1* G 1 -10.681 -1.907 12.234 1.00 0.00 H \n" -"ATOM 18 O2* G 1 -10.882 -0.790 10.004 1.00 0.00 O \n" -"ATOM 19 HO2' G 1 -10.518 -0.117 10.584 1.00 0.00 H \n" -"ATOM 20 N9 G 1 -9.145 -3.330 11.991 1.00 0.00 N \n" -"ATOM 21 C8 G 1 -8.821 -4.645 11.961 1.00 0.00 C \n" -"ATOM 22 H8 G 1 -9.493 -5.380 11.521 1.00 0.00 H \n" -"ATOM 23 N7 G 1 -7.681 -4.993 12.454 1.00 0.00 N \n" -"ATOM 24 C5 G 1 -7.164 -3.763 12.870 1.00 0.00 C \n" -"ATOM 25 C4 G 1 -8.055 -2.738 12.589 1.00 0.00 C \n" -"ATOM 26 C6 G 1 -5.926 -3.463 13.497 1.00 0.00 C \n" -"ATOM 27 N3 G 1 -7.907 -1.418 12.839 1.00 0.00 N \n" -"ATOM 28 C2 G 1 -6.740 -1.143 13.429 1.00 0.00 C \n" -"ATOM 29 N1 G 1 -5.797 -2.102 13.744 1.00 0.00 N \n" -"ATOM 30 H1 G 1 -4.939 -1.812 14.187 1.00 0.00 H \n" -"ATOM 31 N2 G 1 -6.432 0.114 13.748 1.00 0.00 N \n" -"ATOM 32 H21 G 1 -5.549 0.320 14.193 1.00 0.00 H \n" -"ATOM 33 H22 G 1 -7.081 0.860 13.544 1.00 0.00 H \n" -"ATOM 34 O6 G 1 -5.036 -4.227 13.814 1.00 0.00 O \n" -"END"; - - const char* uudPdbString2 = "" -"ATOM 5 C5* G 1 0.000 1.000 0.000 1.00 0.00 C \n" -"ATOM 8 C4* G 1 0.000 0.000 0.000 1.00 0.00 C \n" -"ATOM 10 O4* G 1 1.000 0.000 0.000 1.00 0.00 O \n" -"ATOM 11 C3* G 1 0.000 0.000 1.000 1.00 0.00 C \n" -"END"; - - // std::ifstream inFileStream( "1UUD.pdb"); - std::istringstream inFileStream(uudPdbString); - - PdbStructure pdbStructure(inFileStream); - RNA myRNA("G");// GGC"); - - Compound::AtomTargetLocations atomTargets = myRNA.createAtomTargets(pdbStructure); - - //for (Compound::BondIndex bondIx(0); bondIx < myRNA.getNumBonds(); ++bondIx) - //{ - // myRNA.setBondMobility(BondMobility::Rigid, bondIx); - //} - - // Four steps to a perfect match - myRNA.matchDefaultAtomChirality(atomTargets); - myRNA.matchDefaultBondLengths(atomTargets); - myRNA.matchDefaultBondAngles(atomTargets); - myRNA.matchDefaultDihedralAngles(atomTargets); - myRNA.matchDefaultTopLevelTransform(atomTargets); - - Real residual = myRNA.getTransformAndResidual(atomTargets).residual; - Transform transform = myRNA.getTransformAndResidual(atomTargets).transform; - - cout << "Transform = " << transform << endl; - - cout << "Residual = " << residual << endl; - - myRNA.writeDefaultPdb(std::cout); - - // Make sure H3* atom is within 0.1 Angstrom of target location - Vec3 finalLocH3 = myRNA.getTopLevelTransform() * myRNA.calcDefaultAtomLocationInCompoundFrame("0/H3*"); - Vec3 initialLocH3 = 0.10 * Vec3(-9.675, -3.764, 9.008); - Vec3 errVec = finalLocH3 - initialLocH3; - Real error = std::sqrt( dot(errVec, errVec) ); - - if (error < 0.01) - { - cout << "PASSED" << endl; - - return 0; - } - - else - { - cout << "FAILED" << endl; - cerr << "Initial H3* location was " << initialLocH3 << endl; - cerr << "Final H3* location was " << finalLocH3 << endl; - - return 1; - } -} - --- a/tests/TestLoadPdb.cpp +++ /dev/null @@ -1,220 +0,0 @@ - -#include "molmodel/internal/Compound.h" - -#include "SimTKmolmodel.h" -#include -#include -#include -#include - -using namespace SimTK; -using namespace std; - -void testExactMatch(const std::string& pdbFileName) -{ - time_t startTime = time(NULL); - - // In order for Biotypes to work correctly, - // amber99 parameters must be loaded first? - CompoundSystem system; - SimbodyMatterSubsystem matter(system); - DuMMForceFieldSubsystem forces(system); - forces.loadAmber99Parameters(); - - cout << "Loading Structure from File..." << endl; - std::ifstream pdbInputStream(pdbFileName.c_str()); - assert(pdbInputStream.is_open()); - - // Check speed of various pdb file matching operations - PdbStructure inputPdbStructure(pdbInputStream); - pdbInputStream.close(); - - time_t endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Creating internal coordinate model..." << endl; - RNA rna(inputPdbStructure); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Noting target atom locations..." << endl; - Compound::AtomTargetLocations atomTargets = - rna.createAtomTargets(inputPdbStructure); - cout << atomTargets.size() << " matches found" << endl; - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Matching chirality..." << endl; - rna.matchDefaultAtomChirality(atomTargets, 0.01); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Matching bond angles..." << endl; - rna.matchDefaultBondAngles(atomTargets); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Matching bond lengths..." << endl; - rna.matchDefaultBondLengths(atomTargets); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Matching dihedral angles..." << endl; - rna.matchDefaultDihedralAngles(atomTargets, Compound::DistortPlanarBonds); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Matching top level transform..." << endl; - rna.matchDefaultTopLevelTransform(atomTargets); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Creating simbody system..." << endl; - system.adoptCompound(rna); - rna.assignBiotypes(); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Creating multibody model..." << endl; - system.modelCompounds(); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Writing PDB file" << endl; - PdbStructure outputPdbStructure(rna); - std::ofstream pdbOutputStream("test1.pdb"); - outputPdbStructure.write(pdbOutputStream); - pdbOutputStream.close(); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Done." << endl; -} - -void testIdealizedMatch(const std::string& pdbFileName) -{ - time_t startTime = time(NULL); - - // In order for Biotypes to work correctly, - // amber99 parameters must be loaded first? - CompoundSystem system; - SimbodyMatterSubsystem matter(system); - DuMMForceFieldSubsystem forces(system); - forces.loadAmber99Parameters(); - - cout << "Loading Structure from File..." << endl; - std::ifstream pdbInputStream(pdbFileName.c_str()); - assert(pdbInputStream.is_open()); - - // Check speed of various pdb file matching operations - PdbStructure inputPdbStructure(pdbInputStream); - pdbInputStream.close(); - - time_t endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Creating internal coordinate model..." << endl; - RNA rna(inputPdbStructure); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Noting target atom locations..." << endl; - Compound::AtomTargetLocations atomTargets = - rna.createAtomTargets(inputPdbStructure); - cout << atomTargets.size() << " matches found" << endl; - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Matching chirality..." << endl; - rna.matchDefaultAtomChirality(atomTargets, 90*Deg2Rad); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Matching dihedral angles..." << endl; - rna.matchDefaultDihedralAngles(atomTargets, Compound::KeepPlanarBonds); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Matching top level transform..." << endl; - rna.matchDefaultTopLevelTransform(atomTargets); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Final fitting..." << endl; - rna.fitDefaultConfiguration(atomTargets, 0.005); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Creating simbody system..." << endl; - system.adoptCompound(rna); - rna.assignBiotypes(); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Creating multibody model..." << endl; - system.modelCompounds(); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Writing PDB file" << endl; - PdbStructure outputPdbStructure(rna); - std::ofstream pdbOutputStream("test2.pdb"); - outputPdbStructure.write(pdbOutputStream); - pdbOutputStream.close(); - - endTime = time(NULL); - cout << endTime - startTime << " s elapsed time" << endl; - startTime = endTime; - - cout << "Done." << endl; -} - -int main(int argc, char *argv[]) -{ - // PDB file name is required as argument - assert(argc > 1); - - const std::string pdbFileName(argv[1]); - - testExactMatch(pdbFileName); - testIdealizedMatch(pdbFileName); -} - --- a/tests/TestPdb.cpp +++ /dev/null @@ -1,195 +0,0 @@ -/* -------------------------------------------------------------------------- * - * SimTK Core: SimTK Simbody(tm) * - * -------------------------------------------------------------------------- * - * This is part of the SimTK Core biosimulation toolkit originating from * - * Simbios, the NIH National Center for Physics-Based Simulation of * - * Biological Structures at Stanford, funded under the NIH Roadmap for * - * Medical Research, grant U54 GM072970. See https://simtk.org. * - * * - * Portions copyright (c) 2006-7 Stanford University and the Authors. * - * Authors: Christopher Bruns * - * Contributors: * - * * - * Permission is hereby granted, free of charge, to any person obtaining a * - * copy of this software and associated documentation files (the "Software"), * - * to deal in the Software without restriction, including without limitation * - * the rights to use, copy, modify, merge, publish, distribute, sublicense, * - * and/or sell copies of the Software, and to permit persons to whom the * - * Software is furnished to do so, subject to the following conditions: * - * * - * The above copyright notice and this permission notice shall be included in * - * all copies or substantial portions of the Software. * - * * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * - * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * - * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * - * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, * - * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR * - * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE * - * USE OR OTHER DEALINGS IN THE SOFTWARE. * - * -------------------------------------------------------------------------- */ - -#include "SimTKmolmodel.h" -#include "molmodel/internal/Pdb.h" - -#include - -using namespace SimTK; -using namespace std; - -#define ASSERT(cond) {SimTK_ASSERT_ALWAYS(cond, "Assertion failed");} - -void testInputMatchesOutput() { - // A very small protein - String inputPdb = ""\ - "ATOM 1 N HIS 19 28.165 29.227 23.618 1.00 91.78 N\n" - "ATOM 2 CA HIS 19 27.004 29.173 22.731 1.00 91.74 C\n" - "ATOM 3 C HIS 19 26.321 27.818 22.666 1.00 81.05 C\n" - "ATOM 4 O HIS 19 25.105 27.739 22.805 1.00 89.56 O\n" - "ATOM 5 CB HIS 19 27.248 29.629 21.270 1.00 98.21 C\n" - "ATOM 6 CG HIS 19 27.954 30.936 21.082 1.00100.00 C\n" - "ATOM 7 ND1 HIS 19 28.852 31.112 20.015 1.00100.00 N\n" - "ATOM 8 CD2 HIS 19 27.882 32.111 21.798 1.00100.00 C\n" - "ATOM 9 CE1 HIS 19 29.310 32.368 20.116 1.00100.00 C\n" - "ATOM 10 NE2 HIS 19 28.753 32.997 21.176 1.00100.00 N\n" - "ATOM 11 N SER 20 27.057 26.778 22.303 1.00 58.48 N\n" - "ATOM 12 CA SER 20 26.412 25.501 22.190 1.00 45.35 C\n" - "ATOM 13 C SER 20 26.089 24.909 23.558 1.00 30.40 C\n" - "ATOM 14 O SER 20 26.808 25.067 24.542 1.00 29.60 O\n" - "ATOM 15 CB SER 20 27.206 24.513 21.344 1.00 49.69 C\n" - "ATOM 16 OG SER 20 26.466 23.310 21.049 1.00 22.13 O\n" - "TER\n" - "END\n"; - - istringstream inStream(inputPdb); - PdbStructure pdbStructure(inStream); - - ostringstream outStream; - pdbStructure.write(outStream); - - ASSERT(inputPdb == outStream.str()); -} - - -// SimTK core bug 653 -// Avoid distorting input geometry -void testMatchDefaultBreaksPlanarity() { - // A very small RNA - String pdbString = ""\ - "ATOM 1 O5' A 1 0.000 0.000 0.000 1.00 0.00 O\n" - "ATOM 2 C5' A 1 -0.731 1.221 0.000 1.00 0.00 C\n" - "ATOM 3 H5' A 1 -0.380 1.863 -0.837 1.00 0.00 H\n" - "ATOM 4 H5'' A 1 -0.570 1.749 0.965 1.00 0.00 H\n" - "ATOM 5 C4' A 1 -2.202 0.927 -0.173 1.00 0.00 C\n" - "ATOM 6 H4' A 1 -2.775 1.880 -0.152 1.00 0.00 H\n" - "ATOM 7 O4' A 1 -2.416 0.262 -1.447 1.00 0.00 O\n" - "ATOM 8 C3' A 1 -2.649 -0.071 0.888 1.00 0.00 C\n" - "ATOM 9 H3' A 1 -1.861 -0.843 1.026 1.00 0.00 H\n" - "ATOM 10 O3' A 1 -2.867 0.610 2.118 1.00 0.00 O\n" - "ATOM 11 C2' A 1 -3.982 -0.569 0.339 1.00 0.00 C\n" - "ATOM 12 H2' A 1 -4.164 -1.607 0.693 1.00 0.00 H\n" - "ATOM 13 C1' A 1 -3.707 -0.583 -1.164 1.00 0.00 C\n" - "ATOM 14 H1' A 1 -4.620 -0.263 -1.711 1.00 0.00 H\n" - "ATOM 15 O2' A 1 -4.991 0.378 0.624 1.00 0.00 O\n" - "ATOM 16 HO2' A 1 -5.049 0.485 1.576 1.00 0.00 H\n" - "ATOM 17 N9 A 1 -3.416 -1.965 -1.639 1.00 0.00 N\n" - "ATOM 18 C8 A 1 -2.272 -2.678 -1.387 1.00 0.00 C\n" - "ATOM 19 H8 A 1 -1.653 -2.510 -0.518 1.00 0.00 H\n" - "ATOM 20 N7 A 1 -2.001 -3.578 -2.293 1.00 0.00 N\n" - "ATOM 21 C5 A 1 -3.039 -3.450 -3.209 1.00 0.00 C\n" - "ATOM 22 C4 A 1 -3.909 -2.469 -2.818 1.00 0.00 C\n" - "ATOM 23 C6 A 1 -3.333 -4.123 -4.406 1.00 0.00 C\n" - "ATOM 24 N3 A 1 -5.033 -2.063 -3.454 1.00 0.00 N\n" - "ATOM 25 C2 A 1 -5.370 -2.657 -4.589 1.00 0.00 C\n" - "ATOM 26 N1 A 1 -4.448 -3.699 -5.015 1.00 0.00 N\n" - "ATOM 1 H2 A 1 -6.262 -2.352 -5.114 1.00 0.00 H\n" - "ATOM 28 N6 A 1 -2.564 -5.104 -4.896 1.00 0.00 N\n" - "ATOM 29 H61 A 1 -1.732 -5.389 -4.399 1.00 0.00 H\n" - "ATOM 30 H62 A 1 -2.817 -5.558 -5.761 1.00 0.00 H\n" - "ATOM 31 P A 1 1.610 0.000 0.000 1.00 0.00 P\n" - "ATOM 32 OP1 A 1 2.103 0.698 1.208 1.00 0.00 O\n" - "ATOM 33 OP2 A 1 2.103 0.698 -1.208 1.00 0.00 O\n" - "ATOM 34 OP3 A 1 2.103 -1.395 0.000 1.00 0.00 O\n"; - istringstream pdbStream(pdbString); - PdbStructure pdbStructure(pdbStream); - - // 1) First reproduce undesired behavior - RNA mol1("A"); - Compound::AtomTargetLocations atomTargets = mol1.createAtomTargets(pdbStructure); - - mol1.matchDefaultAtomChirality(atomTargets); - mol1.matchDefaultBondLengths(atomTargets); - mol1.matchDefaultBondAngles(atomTargets); - mol1.matchDefaultDihedralAngles(atomTargets, Compound::DistortPlanarBonds); - mol1.matchDefaultTopLevelTransform(atomTargets); - - // Most distorted part in problem report is distance between atoms N1 and C2 - const ResidueInfo& res1 = mol1.getResidue( ResidueInfo::Index(0) ); - Vec3 atomN1Pos = mol1.calcDefaultAtomFrameInCompoundFrame(res1.getAtomIndex("N1")).p(); - Vec3 atomC2Pos = mol1.calcDefaultAtomFrameInCompoundFrame(res1.getAtomIndex("C2")).p(); - Real bondLength = (atomN1Pos - atomC2Pos).norm(); - - ASSERT(bondLength > 0.20); // distorted - - - // 2) Repair with extra parameter on matchDefaultAtomChirality - - RNA mol2("A"); - atomTargets = mol2.createAtomTargets(pdbStructure); - - mol2.matchDefaultAtomChirality(atomTargets, 0.20); - mol2.matchDefaultBondLengths(atomTargets); - mol2.matchDefaultBondAngles(atomTargets); - mol2.matchDefaultDihedralAngles(atomTargets, Compound::DistortPlanarBonds); - mol2.matchDefaultTopLevelTransform(atomTargets); - - // Most distorted part in problem report is distance between atoms N1 and C2 - const ResidueInfo& res2 = mol2.getResidue( ResidueInfo::Index(0) ); - atomN1Pos = mol2.calcDefaultAtomFrameInCompoundFrame(res2.getAtomIndex("N1")).p(); - atomC2Pos = mol2.calcDefaultAtomFrameInCompoundFrame(res2.getAtomIndex("C2")).p(); - bondLength = (atomN1Pos - atomC2Pos).norm(); - - cout << "bondLength = " << bondLength << endl; - mol2.writeDefaultPdb(cout); - - ASSERT(bondLength < 0.20); // not distorted - - - // 3) Repair by not setting torsion of bonded planar atoms - - RNA mol3("A"); - atomTargets = mol3.createAtomTargets(pdbStructure); - - mol3.matchDefaultAtomChirality(atomTargets); - mol3.matchDefaultBondLengths(atomTargets); - mol3.matchDefaultBondAngles(atomTargets); - mol3.matchDefaultDihedralAngles(atomTargets, Compound::KeepPlanarBonds); - mol3.matchDefaultTopLevelTransform(atomTargets); - - // Most distorted part in problem report is distance between atoms N1 and C2 - const ResidueInfo& res3 = mol3.getResidue( ResidueInfo::Index(0) ); - atomN1Pos = mol3.calcDefaultAtomFrameInCompoundFrame(res3.getAtomIndex("N1")).p(); - atomC2Pos = mol3.calcDefaultAtomFrameInCompoundFrame(res3.getAtomIndex("C2")).p(); - bondLength = (atomN1Pos - atomC2Pos).norm(); - - ASSERT(bondLength < 0.20); // not distorted -} - -int main() { - try { - testInputMatchesOutput(); - testMatchDefaultBreaksPlanarity(); - } - catch (const std::exception& e) { - printf("EXCEPTION THROWN: %s\n", e.what()); - return 1; - } - catch (...) { - printf("UNKNOWN EXCEPTION THROWN\n"); - return 1; - } - - cout << "Done" << endl; - return 0; -} -