Only in src/gleap_new/freelib/readline: libreadline.a Only in src/gleap_new/leapsrc: sleap diff -aur src/gleap/mortsrc/ambfmt/prep.cpp src/gleap_new/mortsrc/ambfmt/prep.cpp --- src/gleap/mortsrc/ambfmt/prep.cpp 2008-06-23 19:43:20.000000000 -0500 +++ src/gleap_new/mortsrc/ambfmt/prep.cpp 2009-03-23 12:20:07.000000000 -0500 @@ -20,7 +20,7 @@ return 10000; } - throw logic_error("Error: unknown chain identifier " + chain); + throw std::runtime_error("Error: unknown chain identifier " + chain); } void read_amber_prep(istream& is, molecule_t& mol) @@ -43,10 +43,10 @@ // the first column, no supported any more. string use_first_column; is >> use_first_column; - if(use_first_column.find("CORR") == string::npos) - { - throw logic_error("Error: prep format, sort by first column no longer supported"); - } + //if(use_first_column.find("CORR") == string::npos) + //{ + // throw std::runtime_error("Error: prep format, sort by first column no longer supported"); + //} is.ignore(MAX_LINE_WIDTH, '\n'); // cutoff: if > 0.0, make bond according to distance @@ -58,6 +58,8 @@ vector crds; vector nunfills; + int head = 0; + int tail = 0; char line[MAX_LINE_WIDTH]; while( is.getline(line, MAX_LINE_WIDTH) && !empty(line) ) { @@ -128,8 +130,18 @@ create_bond(atom, atoms.back() ).set_i(ORDER, 1); nunfills.back() -= 1; } + + if( chain=="M" ) + { + // head is the first M atom, tail is the last M atom + if( head==0 ) + head = atom.absid()+1; + tail = atom.absid() + 1; + } } + + if( nunfills.size()>0 && nunfills.back()==0 ) { crds.pop_back(); @@ -148,6 +160,11 @@ } + m.set_i( HEAD, head ); + m.set_i( TAIL, tail ); + r.set_i( HEAD, head ); + r.set_i( TAIL, tail ); + string keyw; while(is >> keyw && keyw != "DONE") { Only in src/gleap_new/mortsrc/common: libcommon.a Only in src/gleap_new/mortsrc: libmort.a Only in src/gleap_new/mortsrc/nablib: lex.mm_options.c diff -aur src/gleap/mortsrc/object/resdfun.cpp src/gleap_new/mortsrc/object/resdfun.cpp --- src/gleap/mortsrc/object/resdfun.cpp 2008-07-01 13:21:05.000000000 -0500 +++ src/gleap_new/mortsrc/object/resdfun.cpp 2009-03-23 12:55:11.000000000 -0500 @@ -1,3 +1,4 @@ +#include #include #include #include @@ -67,6 +68,9 @@ build_atom(resd, *a); } + std::cout << "mol.natom(): " << mol.natom() << std::endl; + std::cout << "sub.natom(): " << sub.natom() << std::endl; + moiter_t b = sub.bond_begin(); for( ; b != sub.bond_end(); ++b ) { diff -aur src/gleap/plugins/geom.cpp src/gleap_new/plugins/geom.cpp --- src/gleap/plugins/geom.cpp 2008-06-23 19:43:21.000000000 -0500 +++ src/gleap_new/plugins/geom.cpp 2009-03-23 12:32:24.000000000 -0500 @@ -1,11 +1,24 @@ +#include #include - - +#include +#include "impose.hpp" namespace amber { using namespace mort; + void rotate( moref_t& atom, const numvec& orig, const numvec& axis, double angl ) + { + assert( atom.cmpid()==ATOM ); + + numvec pos = atom.get_v(POSITION); + pos -= orig; + mort::rotate( pos, axis, angl ); + pos += orig; + atom.set_v(POSITION, pos); + } + + class geom_finder : public finder_t { public: @@ -117,7 +130,7 @@ */ } - bool impose_dist( const atomvec_t& atoms, double dist ) + void impose_dist( const atomvec_t& atoms, double dist ) { assert( atoms.size() == 2 ); @@ -150,7 +163,9 @@ atomvec_t::iterator ae = target->end(); for( ; ai != ae; ++ai ) { - translate_atom( *ai, offset ); + numvec pos = ai->get_v(POSITION); + pos += offset; + ai->set_v(POSITION, pos ); } /* @@ -163,52 +178,86 @@ } std::cout << std::endl; */ - return true; } - bool impose_angl( const atomvec_t& atoms, double angl ) + + void impose_angl( const numvec& v0, const numvec& v1, const numvec& v2, double angl, atomvec_t& set ) { - assert( atoms.size() == 3 ); + numvec axis = cross( v0-v1, v2-v1 ); + numvec orig = v1; - atomvec_t set_0, set_2; + double curt = mort::angle( v0, v1, v2 ); + double offset = angl - curt; + + atomvec_t::iterator ai = set.begin(); + atomvec_t::iterator ae = set.end(); + for( ; ai != ae; ++ai ) + { + rotate( *ai, orig, axis, offset ); + } - // std::cout << "imposing angl between: " << atoms[0].get_i(ID) << " " << atoms[1].get_i(ID) << " " << atoms[2].get_i(ID) << std::endl; + /* + std::cout << "rotating: " << std::endl; + for( int i=0; i < target->size(); ++i) + { + std::cout << target->at(i).get_i(ID) << " to "; + numvec pos = target->at(i).get_v(POSITION); + std::cout << pos[0] << " " << pos[1] << " " << pos[2] << std::endl; + } + std::cout << std::endl; + */ + } - find_related( atoms[0], atoms[1], set_0 ); - find_related( atoms[2], atoms[1], set_2 ); - numvec v0 = atoms[0].get_v(POSITION); - numvec v1 = atoms[1].get_v(POSITION); - numvec v2 = atoms[2].get_v(POSITION); + void impose_angl( const moref_t& a0, const moref_t& a1, const moref_t& a2, double angl, atomvec_t& set ) + { + numvec v0 = a0.get_v(POSITION); + numvec v1 = a1.get_v(POSITION); + numvec v2 = a2.get_v(POSITION); + impose_angl( v0, v1, v2, angl, set ); + } - numvec axis = cross( v0 - v1, v2 - v1 ); - numvec orig = v1; - double curt = angle( v0, v1, v2 ); + void impose_angl( const atomvec_t& atoms, double angl ) + { + assert( atoms.size()==3 ); + + atomvec_t set_0, set_2; + // std::cout << "imposing angl between: "; + // std::cout << atoms[0].get_i(ID) << " "; + // std::cout << atoms[1].get_i(ID) << " "; + // std::cout << atoms[2].get_i(ID) << std::endl; - atomvec_t* target; - double offset; + find_related( atoms[0], atoms[1], set_0 ); + find_related( atoms[2], atoms[1], set_2 ); + + numvec rot(4); if( set_0.size() <= set_2.size() ) { - target = &set_0; - offset = curt - angl; + impose_angl( atoms[2], atoms[1], atoms[0], angl, set_0 ); } else { - target = &set_2; - offset = angl - curt; + impose_angl( atoms[0], atoms[1], atoms[2], angl, set_2 ); } + } + + void impose_tors( const numvec& v0, const numvec& v1, const numvec& v2, const numvec& v3, double tors, atomvec_t& set ) + { + numvec axis = normalize_copy( v1 - v2 ); + numvec orig = v1; + double curt = mort::torsion( v0, v1, v2, v3 ); + double offset = curt - tors; - atomvec_t::iterator ai = target->begin(); - atomvec_t::iterator ae = target->end(); + atomvec_t::iterator ai = set.begin(); + atomvec_t::iterator ae = set.end(); for( ; ai != ae; ++ai ) { - rotate_atom( *ai, orig, axis, offset ); + rotate( *ai, orig, axis, offset ); } - /* std::cout << "rotating: " << std::endl; for( int i=0; i < target->size(); ++i) @@ -219,19 +268,30 @@ } std::cout << std::endl; */ - - return true; } + void impose_tors( const moref_t& a0, const moref_t& a1, const moref_t& a2, const moref_t& a3, double tors, atomvec_t& set) + { + + std::cout << "imposing tors between: "; + std::cout << a0.get_s(NAME) << "-"; + std::cout << a1.get_s(NAME) << "-"; + std::cout << a2.get_s(NAME) << "-"; + std::cout << a3.get_s(NAME) << std::endl; + + numvec v0 = a0.get_v(POSITION); + numvec v1 = a1.get_v(POSITION); + numvec v2 = a2.get_v(POSITION); + numvec v3 = a3.get_v(POSITION); + + impose_tors( v0, v1, v2, v3, tors, set ); + } void impose_tors( const atomvec_t& atoms, double tors, bool middle_start ) { assert( atoms.size() == 4 ); - // std::cout << "imposing tors between: " << atoms[0].get_i(ID) << " " << atoms[1].get_i(ID) << " " << atoms[2].get_i(ID) << " " << atoms[3].get_i(ID) << std::endl; - - atomvec_t set_1, set_2; if( middle_start ) @@ -245,49 +305,14 @@ find_related( atoms[3], atoms[2], set_2 ); } - numvec v0 = atoms[0].get_v(POSITION); - numvec v1 = atoms[1].get_v(POSITION); - numvec v2 = atoms[2].get_v(POSITION); - numvec v3 = atoms[3].get_v(POSITION); - - numvec axis = normalize_copy( v2 - v1 ); - numvec orig = v1; - double curt = torsion( v0, v1, v2, v3 ); - - atomvec_t* target; - double offset; - - if( set_1.size() <= set_2.size() ) + if( set_1.size() < set_2.size() ) { - target = &set_1; - offset = curt - tors; + impose_tors( atoms[3], atoms[2], atoms[1], atoms[0], tors, set_1 ); } else { - target = &set_2; - offset = tors - curt; + impose_tors( atoms[0], atoms[1], atoms[2], atoms[3], tors, set_2 ); } - - atomvec_t::iterator ai = target->begin(); - atomvec_t::iterator ae = target->end(); - for( ; ai != ae; ++ai ) - { - rotate_atom( *ai, orig, axis, offset ); - } - - - /* - std::cout << "rotating: " << std::endl; - for( int i=0; i < target->size(); ++i) - { - std::cout << target->at(i).get_i(ID) << " to "; - numvec pos = target->at(i).get_v(POSITION); - std::cout << pos[0] << " " << pos[1] << " " << pos[2] << std::endl; - } - std::cout << std::endl; - */ - - } void impose_dist( atomvec_t& atoms, const vector& inter ) diff -aur src/gleap/plugins/impose.cpp src/gleap_new/plugins/impose.cpp --- src/gleap/plugins/impose.cpp 2007-05-09 14:46:06.000000000 -0500 +++ src/gleap_new/plugins/impose.cpp 2009-03-23 12:34:37.000000000 -0500 @@ -64,7 +64,7 @@ } else { - throw logic_error( "Error: can't understand internal constraint " + m_inters[i][0] ); + throw std::runtime_error( "Error: can't understand internal constraint " + m_inters[i][0] ); } } @@ -91,7 +91,7 @@ void impose_command::undo( ) { - throw logic_error( "Sorry, not implemented yet" ); + throw std::runtime_error( "Sorry, not implemented yet" ); } const char* impose_command::info( ) const @@ -103,7 +103,7 @@ { if( args.size() != 4 ) { - throw logic_error( "Error: wrong number of arguments " ); + throw std::runtime_error( "Error: wrong number of arguments " ); } return shared_ptr< command_i >( new impose_command( args[1], interpret1d( args[2] ), interpret2d( args[3] ) ) ); diff -aur src/gleap/plugins/impose.hpp src/gleap_new/plugins/impose.hpp --- src/gleap/plugins/impose.hpp 2008-06-23 19:43:21.000000000 -0500 +++ src/gleap_new/plugins/impose.hpp 2009-03-23 12:34:33.000000000 -0500 @@ -1,6 +1,6 @@ #ifndef GLEAP_PLUGINS_IMPOSE_HPP #define GLEAP_PLUGINS_IMPOSE_HPP - +#include #include #include @@ -40,8 +40,16 @@ void impose_angl( const atomvec_t& atoms, double angl ); + void impose_angl( const numvec& v0, const numvec& v1, const numvec& v2, double angl, atomvec_t& set ); + + void impose_angl( const moref_t& a0, const moref_t& a1, const moref_t& a2, double angl, atomvec_t& set ); + void impose_tors( const atomvec_t& atoms, double tors, bool middle_start ); + void impose_tors( const numvec& v0, const numvec& v1, const numvec& v2, const numvec& v3, double tors, atomvec_t& set ); + + void impose_tors( const moref_t& a0, const moref_t& a1, const moref_t& a2, const moref_t& a3, double tors, atomvec_t& set ); + void impose_dist( atomvec_t& atoms, const vector< string >& inter ); void impose_angl( atomvec_t& atoms, const vector< string >& inter ); Only in src/gleap_new/plugins: libplugins.a diff -aur src/gleap/plugins/merge.cpp src/gleap_new/plugins/merge.cpp --- src/gleap/plugins/merge.cpp 2008-06-23 19:43:21.000000000 -0500 +++ src/gleap_new/plugins/merge.cpp 2009-03-23 12:53:14.000000000 -0500 @@ -3,6 +3,7 @@ #include #include #include "merge.hpp" +#include "impose.hpp" namespace amber { @@ -24,38 +25,108 @@ return n1==a2 ? n2 : n1; } + + moref_t get_nbr( const moref_t& a ) + { + assert( a.related_atom_number() > 0 ); + moiter_t ai = a.related_atom_begin(); + moiter_t ae = a.related_atom_end(); + for( ; ai != ae; ++ai ) + { + if( ai->get_i(ELEMENT) != HYDROGEN ) + return *ai; + } + return a.related_atoms()[0]; + } - atomvec_t get_merge_pts( moref_t& prev ) + numvec get_ang_ref( const moref_t& a ) { + assert( a.related_atom_number() > 0 ); + + numvec orig = a.get_v(POSITION); + + std::vector refs; + moiter_t ai = a.related_atom_begin(); + moiter_t ae = a.related_atom_end(); + for( ; ai != ae; ++ai ) + { + numvec pos = ai->get_v(POSITION); + refs.push_back( normalize_copy(pos - orig) ); + } + + assert( refs.size()==2 || refs.size()==3 ); + + numvec tmp = rotate_copy( refs[0], refs[1], 120.0 ); + if( refs.size()==2 || dist2(tmp,refs[2]) > 0.25) + { + return orig+tmp; + } + + return orig+rotate_copy(refs[0], refs[1], 240.0); + + } + + numvec get_tor_ref( const moref_t& a ) + { + assert( a.related_atom_number() > 0 ); + + std::vector refs; + moiter_t ai = a.related_atom_begin(); + moiter_t ae = a.related_atom_end(); + for( ; ai != ae; ++ai ) + { + numvec pos = ai->get_v(POSITION); + if( ai->get_i(ELEMENT)!=HYDROGEN ) + { + refs.push_back( ai->get_v(POSITION) ); + } + } + + if( refs.size()==0 ) + return a.related_atoms()[0].get_v(POSITION); + + if( refs.size()==1 ) + return refs[0]; + + numvec orig = a.get_v(POSITION); + numvec v0 = normalize_copy( refs[0] - orig ); + numvec v1 = normalize_copy( refs[1] - orig ); + numvec vm = (v0 + v1)*0.5; + return orig + vm; + } + + + + + atomvec_t get_merge_pts( moref_t& prev, moref_t& tail ) + { + atomvec_t atoms; - int tail; - assert( prev.get_i(TAIL, tail) ); - moref_t atm_1( prev.getmol(), ATOM, tail-1 ); - atoms.push_back(atm_1); + atoms.push_back(tail); - if( atm_1.related_atom_number()==0 ) + if( tail.related_atom_number()==0 ) return atoms; - if( atm_1.related_atom_number()==1 ) + if( tail.related_atom_number()==1 ) { - moref_t atm_2 = atm_1.related_atoms()[0]; + moref_t atm_2 = tail.related_atoms()[0]; atoms.push_back( atm_2 ); if( atm_2.related_atom_number()==1 ) return atoms; - moref_t atm_3 = get_other_nbr( atm_2, atm_1 ); + moref_t atm_3 = get_other_nbr( atm_2, tail ); atoms.push_back( atm_3 ); return atoms; } - assert( atm_1.related_atom_number()>=2 ); + assert( tail.related_atom_number()>=2 ); - moref_t atm_2 = atm_1.related_atoms()[0]; - moref_t atm_3 = atm_1.related_atoms()[1]; + moref_t atm_2 = tail.related_atoms()[0]; + moref_t atm_3 = tail.related_atoms()[1]; if( atm_2.get_s(NAME)=="CA" ) { atoms.push_back( atm_2 ); @@ -77,9 +148,11 @@ static const double THETA = DEGRAD*M_PI/180.0; vector crds; - for( int i=0; i < (int)atoms.size(); ++i ) + for( unsigned int i=0; i < atoms.size(); ++i ) { crds.push_back( atoms[i].get_v(POSITION) ); + + //std::cout << "name,pos: " << atoms[i].get_s(NAME) << " " << crds[i][0] << " " << crds[i][1] << " " << crds[i][2] << std::endl; } assert( crds.size()>0 ); @@ -98,28 +171,36 @@ - void place( moref_t& prev, moref_t& curt ) + void place( moref_t& prev, moref_t& curt, moref_t& tail, moref_t& head ) { assert( prev.cmpid()==RESD && curt.cmpid()==RESD ); - numvec pos = get_merge_pos( get_merge_pts(prev) ); - - sparm_comper1 name_is_n (NAME, "N" ); - - int head; - assert( curt.get_i(HEAD, head) ); - moref_t ahead( prev.getmol(), ATOM, head-1 ); + atomvec_t atms; + for( int i=curt.relid(); i < curt.getmol().nresd(); ++i ) + { + atms.push_atom( curt.getmol().resds()[i] ); + } - numvec offset= pos - ahead.get_v(POSITION); + numvec hpos = get_merge_pos( get_merge_pts(prev,tail) ); - moiter_t atom = curt.related_atom_begin(); - for( ; atom != curt.related_atom_end(); ++atom ) + numvec offset = hpos - head.get_v(POSITION); + // translate( curt, offset ); + for(unsigned int i=0; i < atms.size(); ++i ) { - numvec pos = atom->get_v(POSITION); - pos += offset; - atom->set_v(POSITION, pos); + numvec v = atms[i].get_v(POSITION); + atms[i].set_v(POSITION, v+offset); } + + + numvec tref = get_nbr( tail ).get_v(POSITION); + numvec tpos = tail.get_v(POSITION); + + numvec ref1 = get_ang_ref( head ); + impose_angl( tpos, hpos, ref1, 0.0, atms ); + + numvec ref2 = get_tor_ref( head ); + impose_tors( tref, tpos, hpos, ref2, 180.0, atms ); } @@ -145,38 +226,53 @@ if( args.size() == 0 ) { - throw logic_error( "Error : can not understand list: " + m_list ); + throw std::runtime_error( "Error : can not understand list: " + m_list ); } - molecule_ptr pmol( new molecule_t( ) ); + molecule_ptr pmol( new molecule_t() ); - molecule_ptr pres = content().get_mol( args[1] ); - assert( pres != NULL ); + moref_t prev( *pmol, RESD, -1); - moref_t prev = merge( *pmol, *pres ); + int head=0, tail=0; + for( int i = 1; i < (int)args.size() - 1; i++ ) + { + std::cout << "Sequence: " << args[i] << std::endl; + molecule_ptr pres = content().get_mol( args[i] ); - leaplog_t::putline( "Sequence: " + args[1] ); + int size = pmol->natom(); + pmol->get_i( TAIL, tail ); + + moref_t resd = merge( *pmol, *pres, -1 ); + std::cout << "pres.natom: " << pres->natom() << std::endl; + std::cout << "pmol.natom: " << pmol->natom() << std::endl; - for( int i = 2; i < (int)args.size() - 1; i++ ) - { - pres = content().get_mol( args[i] ); - assert( pres != NULL ); - moref_t resd = merge( *pmol, *pres ); - leaplog_t::putline( "Sequence: " + args[i] ); + pres->get_i( HEAD, head ); - int head, tail; - if( m_action == "sequence" && - prev.get_i(TAIL, tail) && tail != 0 && - resd.get_i(HEAD, head) && head != 0 ) + if( m_action == "sequence" && i > 1) { - place( prev, resd ); - moref_t tail_atom( *pmol, ATOM, tail - 1 ); - moref_t head_atom( *pmol, ATOM, head - 1); - moref_t bond = create_bond( tail_atom, head_atom ); - bond.set_i(ORDER, 1); + if( tail==0 ) + { + std::cout << " Tail atom was not set" << std::endl; + } + else if( head==0 ) + { + std::cout << " Head atom was not set" << std::endl; + } + else + { + moref_t tail_atom( *pmol, ATOM, tail-1); + moref_t head_atom( *pmol, ATOM, size+(head-1) ); + std::cout << " tail: " << tail_atom.get_i(ID) << " " << tail_atom.get_s(NAME) << std::endl; + std::cout << " head: " << head_atom.get_i(ID) << " " << head_atom.get_s(NAME) << std::endl; + + place( prev, resd, tail_atom, head_atom ); + moref_t bond = create_bond( tail_atom, head_atom ); + bond.set_i(ORDER, 1); + } } + pres->get_i( TAIL, tail ); prev = resd; } @@ -188,7 +284,7 @@ void merge_command::undo( ) { - throw logic_error( "sorry, not implemented yet." ); + throw std::runtime_error( "sorry, not implemented yet." ); } const char* merge_command::info( ) const @@ -200,7 +296,7 @@ { if( args.size() != 3 ) { - throw logic_error( "Error: wrong number of arguments" ); + throw std::runtime_error( "Error: wrong number of arguments" ); } return shared_ptr< command_i >( new merge_command( args[0], args[1], args[2] ) );