********>Master Bugfix File - AmberTools 1.0 bugfixes 1 to 7. Programs: AmberTools Description: This is the master bugfix file for AmberTools-1.0. It contains all of the bugfixes released against version 1.0 in a single patch script. Usage: Save this file in your $AMBERHOME directory and then apply this patch file to your distribution as follows: cd $AMBERHOME patch -p0 -N -r patch_rejects bugfix 1. Author: Dave Case Date: 7 April 2008 Programs: NAB Severity: major for MPI Description: The nblist.c function has some undefined symbols when MPI is requested. Fix: Apply the following fix to amber10/src/nab/nblist.c Index: nblist.c =================================================================== RCS file: /thr/loyd/case/cvsroot/amber10/src/nab/nblist.c,v retrieving revision 1.5 diff -u -r1.5 nblist.c --- src/nab/nblist.c 29 Mar 2008 19:49:02 -0000 1.5 +++ src/nab/nblist.c 7 Apr 2008 21:18:25 -0000 @@ -950,8 +950,8 @@ threadnum = mycol; numthreads = npcol; } else if (derivs > 0) { - threadnum = mytaskid; - numthreads = numtasks; + threadnum = get_mytaskid(); + numthreads = get_numtasks(); } else { threadnum = 0; numthreads = 1; @@ -963,8 +963,8 @@ threadnum = 0; numthreads = 1; } else { - threadnum = mytaskid; - numthreads = numtasks; + threadnum = get_mytaskid(); + numthreads = get_numtasks(); } #endif --------------------------------------------------------------------------- Workarounds: none. ********> bugfix 2. Author: Dave Case Date: 8 April 2008 Programs: NAB Severity: major for MPI Description: The readparm() routine also doesn't work correctly for MPI Fix: Apply the following fix to amber10/src/nab/prm.c diff -c -r5.20 prm.c *** src/nab/prm.c 17 Mar 2008 20:43:16 -0000 5.20 --- src/nab/prm.c 8 Apr 2008 21:09:43 -0000 *************** *** 819,829 **** #else fscanf(file, " %f", &prm->Pn[i]); #endif ! } ! if( prm->Pn[i] == 0 ){ ! fprintf( stderr, ! "Found an invalid periodicity in the prmtop file: %d\n", i ); ! exit(1); } } --- 819,829 ---- #else fscanf(file, " %f", &prm->Pn[i]); #endif ! if( prm->Pn[i] == 0 ){ ! fprintf( stderr, ! "Found an invalid periodicity in the prmtop file: %d\n", i ); ! exit(1); ! } } } --------------------------------------------------------------------------- Workarounds: none. ********> bugfix 3. Author: Wei Zhang Date: 8 April 2008 Programs: sleap Severity: minor Description: The program won't read a user leaprc file from an arbitray path Fix: Apply the following fix to amber10/src/gleap/mortsrc/guilib/mainwin.cpp *** src/gleap/mortsrc/guilib/mainwin.cpp 30 Aug 2007 18:12:12 -0000 1.10 --- src/gleap/mortsrc/guilib/mainwin.cpp 2 Apr 2008 18:25:11 -0000 1.11 *************** *** 20,25 **** --- 20,33 ---- string find_file( const string& name ) { + std::ifstream is( name.c_str() ); + if( is ) + { + leaplog_t::putline( "using file " + name ); + return name; + } + + string path; assert( mortenv().get_s( "path", path ) ); --------------------------------------------------------------------------- Workarounds: put your leaprc file into one of the standard locations. ********> bugfix 4. Author: Wei Zhang Date: 9 April 2008 Programs: sleap Severity: medium Description: The program doesn't recognize the closeness parameter in solvateBox or solvateOct. Fix: Apply the following fix in amber10/src/gleap. Note that this may make the "solvate" test cases appear to fail. New test cases will be in the next release. *** src/gleap/mortsrc/capbox.hpp 18 Dec 2006 14:54:52 -0000 1.4 --- src/gleap/mortsrc/capbox.hpp 9 Apr 2008 22:48:48 -0000 *************** namespace mort *** 12,18 **** /// \param solvent solvent will be used /// \param center center of the cap /// \param radius radius of the cap ! void solvate_cap( molecule_t& mol, const molecule_t& solvent, const numvec& center, double radius ); /// \brief add a solvation cap to a molecule --- 12,18 ---- /// \param solvent solvent will be used /// \param center center of the cap /// \param radius radius of the cap ! void solvate_cap( molecule_t& mol, const molecule_t& solvent, const numvec& center, double radius, double closeness ); /// \brief add a solvation cap to a molecule *************** namespace mort *** 20,38 **** /// \param solvent solvent will be used /// \param center center of the cap /// \param radius radius of the cap ! void solvate_box( molecule_t& mol, const molecule_t& solvent, double length ); /// \brief solvate a molecule in a box /// \param mol target molecule /// \param solvent solvent will be used /// \param length the box length ! void solvate_oct( molecule_t& mol, const molecule_t& solvent, double length ); /// \brief add a solvation shell to a molecule /// \param mol target molecule /// \param solvent solvent will be used /// \param extent extent of the shell ! void solvate_shell( molecule_t& mol, const molecule_t& solvent, double extent ); /// \brief add an ion to to a molecule /// \param mol target molecule --- 20,38 ---- /// \param solvent solvent will be used /// \param center center of the cap /// \param radius radius of the cap ! void solvate_box( molecule_t& mol, const molecule_t& solvent, double length, double closeness ); /// \brief solvate a molecule in a box /// \param mol target molecule /// \param solvent solvent will be used /// \param length the box length ! void solvate_oct( molecule_t& mol, const molecule_t& solvent, double length, double closeness ); /// \brief add a solvation shell to a molecule /// \param mol target molecule /// \param solvent solvent will be used /// \param extent extent of the shell ! void solvate_shell( molecule_t& mol, const molecule_t& solvent, double extent, double closeness ); /// \brief add an ion to to a molecule /// \param mol target molecule Index: mortsrc/capbox/solvate.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/capbox/solvate.cpp,v retrieving revision 1.16 diff -p -r1.16 solvate.cpp *** src/gleap/mortsrc/capbox/solvate.cpp 21 Jan 2008 22:21:22 -0000 1.16 --- src/gleap/mortsrc/capbox/solvate.cpp 9 Apr 2008 22:48:48 -0000 *************** namespace mort *** 75,98 **** } } - /* - const double* extract_vptr( const molecule_t& mol, const hashid_t& parmid ) - { - assert( mol.get_component(ATOM)->isclean() ); - return mol.get_component(ATOM)->get_vptr(parmid,0); - } - - const double* extract_dptr( const molecule_t& mol, const hashid_t& parmid ) - { - assert( mol.get_component(ATOM)->isclean() ); - return &(mol.get_component(ATOM)->get_d(parmid, 0)); - } - */ - namespace capbox { ! void solvate( molecule_t& mol, const molecule_t& solvent, const numvec& center, const numvec& extent, int shape, double shell_extent ) { numvec box = solvent.get_v(BOX); --- 75,84 ---- } } namespace capbox { ! void solvate( molecule_t& mol, const molecule_t& solvent, const numvec& center, const numvec& extent, int shape, double shell_extent, double closeness ) { numvec box = solvent.get_v(BOX); *************** namespace mort *** 101,137 **** int nz = int( extent[2] / box[2] ) + 1; numvec corner(3); ! corner[0] = center[0] - 0.5*box[0]*nx; ! corner[1] = center[1] - 0.5*box[1]*ny; ! corner[2] = center[2] - 0.5*box[2]*nz; check_shape shaper( shape, center, extent ); ! // std::cout << "nx,ny,nz:"; ! // std::cout << format( "%8d" ) % nx; ! // std::cout << format( "%8d" ) % ny; ! // std::cout << format( "%8d" ) % nz; ! // std::cout << std::endl; int natom = mol.natom(); ! for( int ix=0; ix < nx; ix++ ) { for( int iy=0; iy < ny; iy++ ) { for( int iz=0; iz < nz; iz++ ) { numvec cnt(3); ! cnt[0] = corner[0] + (ix+0.5)*box[0]; ! cnt[1] = corner[1] + (iy+0.5)*box[1]; ! cnt[2] = corner[2] + (iz+0.5)*box[2]; ! ! // std::cout << "place solvent on center: "; ! // std::cout << format( "%8.3f" ) % cnt[0]; ! // std::cout << format( "%8.3f" ) % cnt[1]; ! // std::cout << format( "%8.3f" ) % cnt[2]; ! // std::cout << std::endl; ! check_collision collide( mol, natom, shell_extent, cnt, box ); add_solvent( mol, solvent, cnt, shaper, collide ); } } --- 87,125 ---- int nz = int( extent[2] / box[2] ) + 1; numvec corner(3); ! corner[0] = center[0] + 0.5*box[0]*nx; ! corner[1] = center[1] + 0.5*box[1]*ny; ! corner[2] = center[2] + 0.5*box[2]*nz; check_shape shaper( shape, center, extent ); ! std::cout << "nx,ny,nz:"; ! std::cout << format( "%8d" ) % nx; ! std::cout << format( "%8d" ) % ny; ! std::cout << format( "%8d" ) % nz; ! std::cout << std::endl; int natom = mol.natom(); ! for( int ix=0; ix <= nx; ix++ ) { for( int iy=0; iy < ny; iy++ ) { for( int iz=0; iz < nz; iz++ ) { numvec cnt(3); ! cnt[0] = corner[0] - (ix+0.5)*box[0]; ! cnt[1] = corner[1] - (iy+0.5)*box[1]; ! cnt[2] = corner[2] - (iz+0.5)*box[2]; ! ! /* ! std::cout << "place solvent on center: "; ! std::cout << format( "%8.3f" ) % cnt[0]; ! std::cout << format( "%8.3f" ) % cnt[1]; ! std::cout << format( "%8.3f" ) % cnt[2]; ! std::cout << std::endl; ! */ ! check_collision collide( mol, natom, shell_extent, closeness, cnt, box ); add_solvent( mol, solvent, cnt, shaper, collide ); } } *************** namespace mort *** 247,263 **** return facx + facy + facz < 0.75; } ! check_collision::check_collision( const molecule_t& mol, int natom, double shell_extent, const numvec& center, const numvec& box ) ! : m_shell_extent( shell_extent ) { double max_solute_r = 0.0; moiter_t a = mol.atom_begin(); for( int i=0; i < natom; ++i, ++a ) { ! if( a->get_d(VDWR) > max_solute_r ) { ! max_solute_r = a->get_d(VDWR); } } --- 235,252 ---- return facx + facy + facz < 0.75; } ! check_collision::check_collision( const molecule_t& mol, int natom, double shell_extent, double closeness, const numvec& center, const numvec& box ) ! : m_shell_extent( shell_extent ), ! m_closeness( closeness ) { double max_solute_r = 0.0; moiter_t a = mol.atom_begin(); for( int i=0; i < natom; ++i, ++a ) { ! if( a->get_d(VDWR)*m_closeness > max_solute_r ) { ! max_solute_r = a->get_d(VDWR)*m_closeness; } } *************** namespace mort *** 287,293 **** m_solute_pos.push_back( pos[i3 ] ); m_solute_pos.push_back( pos[i3+1] ); m_solute_pos.push_back( pos[i3+2] ); ! m_solute_vdwr.push_back( a->get_d(VDWR) ); } } --- 276,282 ---- m_solute_pos.push_back( pos[i3 ] ); m_solute_pos.push_back( pos[i3+1] ); m_solute_pos.push_back( pos[i3+2] ); ! m_solute_vdwr.push_back( a->get_d(VDWR)*m_closeness ); } } *************** namespace mort *** 299,305 **** int size = m_solute_vdwr.size(); for( int i=0; i < size; ++i ) { ! double rsum = m_solute_vdwr[i] + rvdw; double rsum2 = rsum*rsum; double dx = m_solute_pos[3*i ] - solventpos[0]; --- 288,294 ---- int size = m_solute_vdwr.size(); for( int i=0; i < size; ++i ) { ! double rsum = m_solute_vdwr[i] + rvdw*m_closeness; double rsum2 = rsum*rsum; double dx = m_solute_pos[3*i ] - solventpos[0]; *************** namespace mort *** 328,379 **** } // namespace capbox - void solvate_cap( molecule_t& mol, const molecule_t& solvent, const numvec& center, double radius ) - { - numvec size( 2*radius, 2*radius, 2*radius ); - capbox::solvate( mol, solvent, center, size, CAP, 0.0 ); - - mol.set_i(SOLUTE, CAP); - mol.set_v(CAP, numvec(center[0], center[1], center[2], radius) ); - } - - void solvate_box( molecule_t& mol, const molecule_t& solvent, double buffer ) - { - set_parm99_vdwr( mol ); - - numvec pmin(3), pmax(3); - minmaxpos( mol, pmin, pmax, true ); - assert( pmin.size()==3 && pmax.size()==3 ); - - numvec center = (pmin + pmax) / 2; - numvec extent = (pmax - pmin); - extent[0] += 2*buffer; - extent[1] += 2*buffer; - extent[2] += 2*buffer; - - std::cout << "box size : "; - std::cout << format( "%10.3f" ) % extent[0]; - std::cout << format( "%10.3f" ) % extent[1]; - std::cout << format( "%10.3f" ) % extent[2]; - std::cout << std::endl; - - - moiter_t a = mol.atom_begin(); - for( ; a != mol.atom_end(); ++a ) - { - numvec p = a->get_v(POSITION); - p -= center; - a->set_v(POSITION, p); - } - - capbox::solvate( mol, solvent, numvec(0.0,0.0,0.0), extent, BOX, 0.0 ); - - mol.set_i(SOLUTE, BOX); - - mol.set_v(BOX, numvec( extent[0], extent[1], extent[2], 90.0 ) ); - } - double ewald_rotate( molecule_t& m ) { int natom = m.natom(); --- 317,323 ---- *************** namespace mort *** 445,458 **** return buffer * tmp / bmax; } ! void solvate_oct( molecule_t& mol, const molecule_t& solvent, double buffer ) { set_parm99_vdwr( mol ); numvec cnt = center( mol ); move( mol, -1.0 * cnt ); - matrix moment = inter_moment( mol ); int nrot; --- 389,447 ---- return buffer * tmp / bmax; } + numvec center_by_vdwr( molecule_t& mol ) + { + numvec pmin(3), pmax(3); + minmaxpos( mol, pmin, pmax, true ); + + numvec center = (pmin + pmax) / 2; + numvec extent = pmax - pmin; + + double* pcrd = get_vptr( mol, ATOM, POSITION ); + for( int i=0; i < mol.natom(); ++i ) + { + pcrd[3*i ] -= center[0]; + pcrd[3*i+1] -= center[1]; + pcrd[3*i+2] -= center[2]; + } + + return extent; + } + + void solvate_cap( molecule_t& mol, const molecule_t& solvent, const numvec& center, double radius, double closeness ) + { + numvec size( 2*radius, 2*radius, 2*radius ); + capbox::solvate( mol, solvent, center, size, CAP, 0.0, closeness ); + + mol.set_i(SOLUTE, CAP); + mol.set_v(CAP, numvec(center[0], center[1], center[2], radius) ); + } + + void solvate_box( molecule_t& mol, const molecule_t& solvent, double buffer, double closeness ) + { + set_parm99_vdwr( mol ); + numvec extent = center_by_vdwr( mol ); + + extent[0] += 2*buffer; + extent[1] += 2*buffer; + extent[2] += 2*buffer; + + capbox::solvate( mol, solvent, numvec(0.0,0.0,0.0), extent, BOX, 0.0, closeness ); ! mol.set_i(SOLUTE, BOX); ! mol.set_v(BOX, numvec( extent[0], extent[1], extent[2], 90.0 ) ); ! ! set_parm99_vdwr( mol ); ! center_by_vdwr( mol ); ! ! } ! ! void solvate_oct( molecule_t& mol, const molecule_t& solvent, double buffer, double closeness ) { set_parm99_vdwr( mol ); numvec cnt = center( mol ); move( mol, -1.0 * cnt ); matrix moment = inter_moment( mol ); int nrot; *************** namespace mort *** 475,485 **** rotate(pos+3*i, rot); } ! numvec pmin(3), pmax(3); ! minmaxpos( mol, pmin, pmax, true ); ! ! cnt = ( pmin + pmax ) / 2.0; ! move( mol, -1.0 * cnt ); double minbuf = min_oct_buf( mol, buffer ); if( buffer < minbuf ) --- 464,470 ---- rotate(pos+3*i, rot); } ! numvec extent = center_by_vdwr( mol ); double minbuf = min_oct_buf( mol, buffer ); if( buffer < minbuf ) *************** namespace mort *** 487,534 **** buffer = minbuf; } ! ! numvec extent = ( pmax - pmin ) + numvec( 2*buffer, 2*buffer, 2*buffer ); double extmax = max( extent ); extent = numvec(extmax, extmax, extmax); ! capbox::solvate( mol, solvent, numvec(0.0, 0.0, 0.0), extent, OCT, 0.0 ); ewald_rotate( mol ); mol.set_i(SOLUTE, BOX); - double boxlen = extmax * sqrt(0.75); mol.set_v(BOX, numvec(boxlen, boxlen, boxlen, 109.471219) ); } ! void solvate_shell( molecule_t& mol, const molecule_t& solvent, double buffer ) { set_parm99_vdwr( mol ); ! numvec pmin(3), pmax(3); ! minmaxpos( mol, pmin, pmax, true ); ! ! numvec center = (pmin + pmax) / 2; ! numvec extent = (pmax - pmin) + numvec( 2*buffer, 2*buffer, 2*buffer ); ! /* ! std::cout << "box size : "; ! std::cout << format( "%10.3f" ) % extent[0]; ! std::cout << format( "%10.3f" ) % extent[1]; ! std::cout << format( "%10.3f" ) % extent[2]; ! std::cout << std::endl; ! ! a = atom_begin(mol); ! for( ; a != atom_end(mol); ++a ) ! { ! numvec p = a.get< position_t >(); ! p -= center; ! a.set< position_t >(p); ! } ! */ ! ! capbox::solvate( mol, solvent, center, extent, BOX, buffer ); } --- 472,504 ---- buffer = minbuf; } ! extent += numvec( 2*buffer, 2*buffer, 2*buffer ); double extmax = max( extent ); extent = numvec(extmax, extmax, extmax); ! capbox::solvate( mol, solvent, numvec(0.0, 0.0, 0.0), extent, OCT, 0.0, closeness ); ! + set_parm99_vdwr( mol ); + center_by_vdwr( mol ); + ewald_rotate( mol ); mol.set_i(SOLUTE, BOX); double boxlen = extmax * sqrt(0.75); mol.set_v(BOX, numvec(boxlen, boxlen, boxlen, 109.471219) ); } ! void solvate_shell( molecule_t& mol, const molecule_t& solvent, double buffer, double closeness ) { set_parm99_vdwr( mol ); + numvec extent = center_by_vdwr( mol ); ! extent += numvec( 2*buffer, 2*buffer, 2*buffer ); ! capbox::solvate( mol, solvent, numvec(0.0, 0.0, 0.0), extent, BOX, buffer, closeness ); + set_parm99_vdwr( mol ); + center_by_vdwr( mol ); } Index: mortsrc/capbox/solvate.hpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/capbox/solvate.hpp,v retrieving revision 1.4 diff -p -r1.4 solvate.hpp *** src/gleap/mortsrc/capbox/solvate.hpp 21 Dec 2006 17:27:15 -0000 1.4 --- src/gleap/mortsrc/capbox/solvate.hpp 9 Apr 2008 22:48:48 -0000 *************** namespace mort *** 47,53 **** struct check_collision { ! check_collision( const molecule_t& mol, int natom, double shell_extent, const numvec& center, const numvec& box ); bool operator()( const double* pos, double rvdw ); --- 47,53 ---- struct check_collision { ! check_collision( const molecule_t& mol, int natom, double shell_extent, double closeness, const numvec& center, const numvec& box ); bool operator()( const double* pos, double rvdw ); *************** namespace mort *** 56,64 **** vector< double > m_solute_vdwr; double m_shell_extent; }; ! void solvate( molecule_t& mol, const molecule_t& solvent, const numvec& center, const numvec& length, int shape, double shell_extent ); void add_solvent( molecule_t& mol, const molecule_t& solvent, const numvec& pos, check_shape shaper, check_collision collide ); --- 56,66 ---- vector< double > m_solute_vdwr; double m_shell_extent; + + double m_closeness; }; ! void solvate( molecule_t& mol, const molecule_t& solvent, const numvec& center, const numvec& length, int shape, double shell_extent, double closeness ); void add_solvent( molecule_t& mol, const molecule_t& solvent, const numvec& pos, check_shape shaper, check_collision collide ); ? plugins/.bond.cpp.swp ? plugins/addmap.cpp.patch ? plugins/makedepend ? plugins/t ? plugins/testa.cpp Index: plugins/solvate.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/plugins/solvate.cpp,v retrieving revision 1.7 diff -p -r1.7 solvate.cpp *** src/gleap/plugins/solvate.cpp 9 May 2007 19:46:06 -0000 1.7 --- src/gleap/plugins/solvate.cpp 9 Apr 2008 22:49:32 -0000 *************** namespace amber *** 11,19 **** } solvate_command::solvate_command( const string& action, const string& solute, const string& solvent, ! const string& center, const string& extent ) : m_action( action ), m_solute( solute ), m_center( center ), m_solvent( solvent ), ! m_extent( extent ) { } --- 11,19 ---- } solvate_command::solvate_command( const string& action, const string& solute, const string& solvent, ! const string& center, const string& extent, double closeness ) : m_action( action ), m_solute( solute ), m_center( center ), m_solvent( solvent ), ! m_extent( extent ), m_closeness( closeness ) { } *************** namespace amber *** 33,39 **** { double length = atof( m_extent.c_str() ); ! solvate_box( *psolute, *psolvent, length ); return true; } --- 33,39 ---- { double length = atof( m_extent.c_str() ); ! solvate_box( *psolute, *psolvent, length, m_closeness ); return true; } *************** namespace amber *** 42,48 **** { double length = atof( m_extent.c_str() ); ! solvate_oct( *psolute, *psolvent, length ); return true; } --- 42,48 ---- { double length = atof( m_extent.c_str() ); ! solvate_oct( *psolute, *psolvent, length, m_closeness ); return true; } *************** namespace amber *** 73,79 **** double radius = atof( m_extent.c_str() ); ! solvate_cap( *psolute, *psolvent, cnt, radius ); return true; } --- 73,79 ---- double radius = atof( m_extent.c_str() ); ! solvate_cap( *psolute, *psolvent, cnt, radius, m_closeness ); return true; } *************** namespace amber *** 82,88 **** double extent = atof( m_extent.c_str() ); ! solvate_shell( *psolute, *psolvent, extent ); return true; } --- 82,88 ---- double extent = atof( m_extent.c_str() ); ! solvate_shell( *psolute, *psolvent, extent, m_closeness ); return true; } *************** namespace amber *** 103,124 **** { if( args[0]=="solvatebox" || args[0]=="solvateshell" || args[0]=="solvateoct" ) { ! if( args.size() != 4 ) { throw logic_error( "Error: wrong number of arguments" ); } ! ! return shared_ptr< command_i >( new solvate_command( args[0], args[1], args[2], "", args[3] ) ); } assert( args[0] == "solvatecap" ); ! if( args.size() != 5 ) { throw logic_error( "Error: wrong number of arguments" ); } ! return shared_ptr< command_i >( new solvate_command( args[0], args[1], args[2], args[3], args[4] ) ); } } // namespace amber --- 103,126 ---- { if( args[0]=="solvatebox" || args[0]=="solvateshell" || args[0]=="solvateoct" ) { ! if( args.size()!=4 && args.size()!=5 ) { throw logic_error( "Error: wrong number of arguments" ); } ! ! double closeness = (args.size()==4) ? 1.0 : atof( args[4].c_str() ); ! return shared_ptr< command_i >( new solvate_command( args[0], args[1], args[2], "", args[3], closeness ) ); } assert( args[0] == "solvatecap" ); ! if( args.size()!=5 && args.size()!=6 ) { throw logic_error( "Error: wrong number of arguments" ); } ! double closeness = (args.size()==5) ? 1.0 : atof( args[5].c_str() ); ! return shared_ptr< command_i >( new solvate_command( args[0], args[1], args[2], args[3], args[4], closeness ) ); } } // namespace amber Index: plugins/solvate.hpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/plugins/solvate.hpp,v retrieving revision 1.2 diff -p -r1.2 solvate.hpp *** src/gleap/plugins/solvate.hpp 9 May 2007 19:46:06 -0000 1.2 --- src/gleap/plugins/solvate.hpp 9 Apr 2008 22:49:32 -0000 *************** namespace amber *** 13,19 **** solvate_command( const string& type ); ! solvate_command( const string& action, const string& solute, const string& solvent, const string& center, const string& extent ); virtual ~solvate_command( ); --- 13,19 ---- solvate_command( const string& type ); ! solvate_command( const string& action, const string& solute, const string& solvent, const string& center, const string& extent, double closeness ); virtual ~solvate_command( ); *************** namespace amber *** 36,41 **** --- 36,43 ---- string m_solvent; string m_extent; + + double m_closeness; }; } // namespace amber --------------------------------------------------------------------------- Workarounds: none ********> bugfix 5. Author: Wei Zhang Date: 24 April 2008 Programs: sleap Severity: Seems to be specific to Mac OSX Description: Segfaults can occur when running sleap Fix: Apply the following fix in amber10/src/gleap. --------------------------------------------------------------------------- =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/leapsrc/Makefile,v retrieving revision 1.17 diff -p -r1.17 Makefile *** src/gleap/leapsrc/Makefile 26 Mar 2008 16:44:26 -0000 1.17 --- src/gleap/leapsrc/Makefile 8 Apr 2008 19:32:19 -0000 *************** INCPATH = -I../mortsrc -I../mortsr *** 21,30 **** ifeq ($(OSTYPE),darwin) LINK = g++ -static-libgcc else ! LINK = g++ -static endif LFLAGS = ! LIBS = ../plugins/libplugins.a ../mortsrc/libmort.a ../freelib/readline/libreadline.a -lm AR = ar cqs RANLIB = QMAKE = /usr/bin/qmake --- 21,30 ---- ifeq ($(OSTYPE),darwin) LINK = g++ -static-libgcc else ! LINK = g++ -static-libgcc endif LFLAGS = ! LIBS = ../freelib/readline/libreadline.a -lm ../plugins/libplugins.a ../mortsrc/libmort.a AR = ar cqs RANLIB = QMAKE = /usr/bin/qmake *************** OBJECTS_DIR = ./ *** 49,60 **** ####### Files ! SOURCES = strbuff.cpp \ ! strmain.cpp \ ! plugins.cpp ! OBJECTS = strbuff.o \ ! strmain.o \ ! plugins.o DESTDIR = TARGET = sleap --- 49,56 ---- ####### Files ! SOURCES = strmain.cpp strbuff.cpp plugins.cpp ! OBJECTS = strmain.o strbuff.o plugins.o DESTDIR = TARGET = sleap Index: leapsrc/plugins.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/leapsrc/plugins.cpp,v retrieving revision 1.11 diff -p -r1.11 plugins.cpp *** src/gleap/leapsrc/plugins.cpp 20 Mar 2008 22:08:41 -0000 1.11 --- src/gleap/leapsrc/plugins.cpp 8 Apr 2008 19:32:19 -0000 *************** *** 32,38 **** --- 32,40 ---- #include #include #include + #include + amber::help_command g_help_command; amber::add_command g_add_command2; amber::addions_command g_addions_command2; amber::addmap_command g_addpdbatommap_command2( "addpdbatommap" ); *************** amber::saveoff_command g_saveoff_command *** 82,85 **** amber::savemol2_command g_savemol2_command2; amber::energy_command g_energy_command2; amber::loadprep_command g_loadprep_command2; - --- 84,86 ---- Index: leapsrc/strmain.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/leapsrc/strmain.cpp,v retrieving revision 1.14 diff -p -r1.14 strmain.cpp *** src/gleap/leapsrc/strmain.cpp 18 Jan 2008 21:54:01 -0000 1.14 --- src/gleap/leapsrc/strmain.cpp 8 Apr 2008 19:32:19 -0000 *************** int main( int argc, char** argv ) *** 59,69 **** - // for windows platform, you need to use at least one symbol defined in shared library, - // to make it be loaded to memory - - #include "../plugins/help.hpp" - - amber::help_command g_help_command; - - --- 59,61 ---- Index: mortsrc/guilib/command.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/guilib/command.cpp,v retrieving revision 1.2 diff -p -r1.2 command.cpp *** src/gleap/mortsrc/guilib/command.cpp 18 Dec 2006 14:48:47 -0000 1.2 --- src/gleap/mortsrc/guilib/command.cpp 8 Apr 2008 19:32:19 -0000 *************** *** 3,19 **** namespace mort { command_i::command_i() { } command_i::command_i( const string& name ) { ! control_t::insert( name, this ); } command_i::~command_i() { } } // namespace amber --- 3,43 ---- namespace mort { + using std::map; + command_i::command_i() { } command_i::command_i( const string& name ) { ! insert( name, this ); } command_i::~command_i() { + } + + map& command_i::dict() + { + static map g_commands; + return g_commands; + } + + void command_i::insert( const string& name, command_i* inst ) + { + dict()[name] = inst; + } + + command_i* command_i::find( const string& name ) + { + map::iterator i = dict().find( name ); + if( i==dict().end() ) + { + throw logic_error( name + ": command not found!" ); + } + + return i->second; } } // namespace amber Index: mortsrc/guilib/command.hpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/guilib/command.hpp,v retrieving revision 1.3 diff -p -r1.3 command.hpp *** src/gleap/mortsrc/guilib/command.hpp 28 Mar 2007 15:32:09 -0000 1.3 --- src/gleap/mortsrc/guilib/command.hpp 8 Apr 2008 19:32:19 -0000 *************** namespace mort *** 30,35 **** --- 30,43 ---- /// making a new command object with given argument. virtual shared_ptr clone( const vector& args ) const = 0; + + static void insert( const string& name, command_i* inst ); + + static command_i* find( const string& name ); + + private: + + static std::map& dict(); }; typedef vector argvec_t; Index: mortsrc/guilib/control.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/guilib/control.cpp,v retrieving revision 1.16 diff -p -r1.16 control.cpp *** src/gleap/mortsrc/guilib/control.cpp 5 Feb 2008 21:42:16 -0000 1.16 --- src/gleap/mortsrc/guilib/control.cpp 8 Apr 2008 19:32:19 -0000 *************** namespace mort *** 131,136 **** --- 131,137 ---- bool control_t::run( const string& command ) { + assert( ! command.empty() ); vector args; *************** namespace mort *** 163,176 **** string name = args[0]; ! map::iterator i = g_commands.find( name ); ! if( i == g_commands.end() ) ! { ! throw logic_error( name + ": command not found!" ); ! } ! command_ptr comm = i->second->clone(args); comm->exec(); --- 164,174 ---- string name = args[0]; ! command_i* pcmd = command_i::find( name ); ! assert( pcmd != NULL ); ! command_ptr comm = pcmd->clone(args); comm->exec(); Index: mortsrc/pdbent/build.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/pdbent/build.cpp,v retrieving revision 1.28 diff -p -r1.28 build.cpp *** src/gleap/mortsrc/pdbent/build.cpp 12 Dec 2007 15:19:41 -0000 1.28 --- src/gleap/mortsrc/pdbent/build.cpp 8 Apr 2008 19:32:19 -0000 *************** namespace mort *** 374,380 **** } else { - //std::cout << "building: " << type << std::endl; pdbent::build_bymdl( *ri, *mdl, nmap, dst, idmap ); } --- 374,379 ---- Index: plugins/add.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/plugins/add.cpp,v retrieving revision 1.4 diff -p -r1.4 add.cpp *** src/gleap/plugins/add.cpp 9 May 2007 19:46:06 -0000 1.4 --- src/gleap/plugins/add.cpp 8 Apr 2008 19:32:19 -0000 *************** namespace amber *** 42,48 **** throw logic_error( "Error: only an atom list is allowed to be added into a unit" ); } ! insert( *unit, *list ); return true; } --- 42,48 ---- throw logic_error( "Error: only an atom list is allowed to be added into a unit" ); } ! mort::insert( *unit, *list ); return true; } *************** namespace amber *** 106,115 **** } amber::add_command g_add_command; - - amber::null_command g_addpdbatommap_command( "addpdbatommap" ); - - amber::null_command g_addpdbresmap_command( "addpdbresmap" ); amber::null_command g_addatomtypes_command( "addatomtypes" ); --- 106,111 ---- Index: plugins/loadpdb.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/plugins/loadpdb.cpp,v retrieving revision 1.15 diff -p -r1.15 loadpdb.cpp *** src/gleap/plugins/loadpdb.cpp 11 Dec 2007 15:46:12 -0000 1.15 --- src/gleap/plugins/loadpdb.cpp 8 Apr 2008 19:32:19 -0000 *************** namespace amber *** 55,60 **** --- 55,62 ---- if( m_action=="loadpdb" ) { + // pff is pointer to force field, it is needed here because for + // inhibitor, we run parmchk to get force field parameters molecule_ptr pff = content().get_mol( "_amber-atom" ); build_bymdl( *pmol, content(), *nmap, pff.get() ); } --------------------------------------------------------------------------- Workarounds: none ********> bugfix 6. Author: Wei Zhang Date: 29 April 2008 Programs: sleap Severity: addions fix for solvated systems Description: Addions doesn't work correctly with solvated systems. Fix: Apply the following fix in amber10/src/gleap. --------------------------------------------------------------------------- Index: mortsrc/capbox/addions.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/capbox/addions.cpp,v retrieving revision 1.10 diff -p -r1.10 addions.cpp *** src/gleap/mortsrc/capbox/addions.cpp 26 Jun 2007 13:44:00 -0000 1.10 --- src/gleap/mortsrc/capbox/addions.cpp 29 Apr 2008 19:20:55 -0000 *************** *** 1,61 **** #include #include "octree.hpp" namespace mort { using namespace capbox; ! void add_ion( molecule_t& mol, const molecule_t& ion, int number, double shell_extent, double resolution ) { numvec pmin(3), pmax(3); ! minmaxpos( mol, pmin, pmax, true ); ! assert( pmin.size()==3 && pmax.size()==3 ); ! double add_extent = ion.atom_begin()->get_d(VDWR); ! add_extent = 1.908; ! pmin[0] = pmin[0] - add_extent - shell_extent; ! pmin[1] = pmin[1] - add_extent - shell_extent; ! pmin[2] = pmin[2] - add_extent - shell_extent; ! ! pmax[0] = pmax[0] + add_extent + shell_extent; ! pmax[1] = pmax[1] + add_extent + shell_extent; ! pmax[2] = pmax[2] + add_extent + shell_extent; ! ! double edge = pow( 2, int( log( max( pmax - pmin ) ) / log( 2 ) + 1.0 ) ); ! ! octree_t tree( pmin, edge ); ! atomvec_t atoms( mol.atoms() ); ! tree.divide( atoms, resolution, checkfun_t( add_extent, shell_extent ) ); ! ! std::pair min = make_pair(0.0, numvec(3)); ! std::pair max = make_pair(0.0, numvec(3)); ! tree.calculate( capbox::elepot_t(mol), min, max ); ! for( int i=0; i < number; ++i ) ! { ! //assert( min.second.size()==3 ); ! //assert( max.second.size()==3 ); numvec pos = charge(ion)>0 ? min.second : max.second; ! /* ! std::cout << "min:"; ! std::cout << min.first << " "; ! std::cout << min.second[0] << " "; ! std::cout << min.second[1] << " "; ! std::cout << min.second[2] << std::endl; ! ! std::cout << "max:"; ! std::cout << max.first << " "; ! std::cout << max.second[0] << " "; ! std::cout << max.second[1] << " "; ! std::cout << max.second[2] << std::endl; ! */ ! ! moref_t rion = merge( mol, ion ); ! moveto(rion, min.second); ! tree.calculate(elepot_t(ion), min, max); ! } } } // namespace mort --- 1,139 ---- + #include #include #include "octree.hpp" + #include "solvate.hpp" namespace mort { + double max_vdwr( molecule_t& m ) + { + double* pr = get_dptr( m, ATOM, VDWR ); + return *std::max_element( pr, pr+m.natom() ); + } + + moref_t find_nearest_solvent( const molecule_t& m, const numvec& pos, int solvent_start=0 ) + { + std::cout << "Info: looking for solvent at: " << pos[0] << " " << pos[1] << " " << pos[2] << std::endl; + + moref_t first( m, ATOM, solvent_start ); + + int solvent_resd_size = first.related_resd_begin()->related_atom_number(); + + const double* pcrd = get_vptr( m, ATOM, POSITION ); + + double mind2 = 1e20; + int minid = -1; + + int natom = m.natom(); + for( int i=solvent_start; i < natom; i+=solvent_resd_size ) + { + double dx = pcrd[3*i ] - pos[0]; + double dy = pcrd[3*i+1] - pos[1]; + double dz = pcrd[3*i+2] - pos[2]; + double d2 = dx*dx + dy*dy + dz*dz; + if( minid < 0 || mind2 > d2 ) + { + mind2 = d2; + minid = i; + } + } + + return moref_t( m, ATOM, minid ); + } + + using namespace capbox; ! void add_ion( molecule_t& mol, const molecule_t& ion, int num_ion, double shell_extent, double resolution ) { + int solute_size = mol.has_i(SOLUTESIZE) ? mol.get_i(SOLUTESIZE) : mol.natom(); + set_amberff_vdwr( mol ); + numvec pmin(3), pmax(3); ! minmaxpos( mol, pmin, pmax, false, solute_size ); ! double max_radius = max_vdwr( mol ); ! double ion_radius = get_amberff_vdwr( ion.atom_begin()->get_s(TYPE) ); ! if( ion.natom() > 1 ) ! { ! std::cout << "Info: polyatomic ion, radius set to " << ion_radius*2.5 << " (2.5 times first atom)" << std::endl; ! ion_radius *= 2.5; ! } ! if( std::abs( charge(ion) ) < 0.1 ) ! { ! throw std::runtime_error( "Error: the ion you chose is in fact netural, thus add ion won't work" ); ! } ! double buffer = ion_radius + max_radius + shell_extent; ! std::cout << "Info: rmax, rion, shell: " << max_radius << " " << ion_radius << " " << shell_extent << std::endl; ! ! pmin -= numvec(buffer, buffer, buffer); ! pmax += numvec(buffer, buffer, buffer); ! std::cout << "Info: enclosing: " << pmin[0] << " " << pmin[1] << " " << pmin[2] << std::endl; ! std::cout << "Info: to: " << pmax[0] << " " << pmax[1] << " " << pmax[2] << std::endl; + + double edge = max( pmax - pmin ); + int ndepth = int( log(edge/resolution)/log(2.0) ) + 1; + octree_t tree( pmin, ndepth, resolution ); + + + double* pcrd = get_vptr( mol, ATOM, POSITION ); + double* prad = get_dptr( mol, ATOM, VDWR ); + + shell_shaper shaper(pcrd, prad, solute_size, ion_radius, shell_extent); + tree.make_shape( shaper ); + + std::pair min = make_pair( 1e20, numvec(3)); + std::pair max = make_pair(-1e20, numvec(3)); + tree.calculate( capbox::elepot_t(mol, solute_size), min, max ); + + vector crds; + + for( int i=0; i < num_ion; ++i ) + { numvec pos = charge(ion)>0 ? min.second : max.second; ! crds.push_back( pos ); ! ! moref_t rion( mol, RESD, -1 ); ! ! if( mol.natom()==solute_size+i ) ! { ! // no solvent, simply add the solvent to the end ! std::cout << "Info: new ion will be placed at: " << pos[0] << " " << pos[1] << " " << pos[2] << " " << std::endl; ! rion = merge( mol, ion ); ! } ! else ! { ! // solvent present, find the nearest solvent residue, use ion to replace it. ! // besides the ion must be before solvent, after the solute ! ! moref_t a = find_nearest_solvent( mol, pos, solute_size ); ! pos = a.get_v(POSITION); ! std::cout << "Info: new ion will be placed at: " << pos[0] << " " << pos[1] << " " << pos[2] << " " << std::endl; ! ! int nsolute_resd = (mol.atom_begin() + solute_size)->related_resd_begin()->relid(); ! rion = merge( mol, ion, nsolute_resd+i ); ! ! mol.remove_resd( *a.related_resd_begin() ); ! //mol.cleanup(); ! } ! ! moveto(rion, pos ); ! if( i < num_ion-1 ) ! { ! // more ions to add: cut a sphere from the tree, make sure no overlapping ! // the radius is 2*ion_radius, to avoid in the future ion ! // with inverse sign will be placed close to this ion. ! outsphere_shaper shaper2(pos, ion_radius*2); ! tree.make_shape( shaper2 ); ! ! min.first = 1e20; ! max.first =-1e20; ! tree.calculate(elepot_t(pos, charge(ion)), min, max); ! } ! } } } // namespace mort Index: mortsrc/capbox/octree.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/capbox/octree.cpp,v retrieving revision 1.5 diff -p -r1.5 octree.cpp *** src/gleap/mortsrc/capbox/octree.cpp 21 Dec 2006 17:27:01 -0000 1.5 --- src/gleap/mortsrc/capbox/octree.cpp 29 Apr 2008 19:20:55 -0000 *************** namespace mort *** 7,182 **** namespace capbox { ! enum node_type { INCLUDED, EXCLUDED, PARTIAL }; ! ! octree_t::octree_t( const numvec& corner, double length ) ! : m_corner(3), m_center(3) { ! m_corner = corner; ! m_length = length; ! m_radius = length * sqrt( 0.75 ); ! ! m_center[0] = m_corner[0] + length / 2.0; ! m_center[1] = m_corner[1] + length / 2.0; ! m_center[2] = m_corner[2] + length / 2.0; } octree_t::~octree_t() { } ! ! int final_check( const numvec& corner, const atomvec_t& atoms, const checkfun_t& checkfun ) { ! atomvec_t::const_iterator ai = atoms.begin(); ! ! int type = EXCLUDED; ! ! for( ; ai != atoms.end(); ++ai ) { ! int pos = checkfun( corner, *ai ); ! ! if( pos == TOTAL_INNER ) ! { ! return EXCLUDED; ! } ! ! if( pos == TOTAL_SHELL ) ! { ! type = INCLUDED; ! } } ! ! return type; ! } ! ! int octree_t::divide( const atomvec_t& atoms, double resolution, const checkfun_t& checkfun ) ! { ! if( m_length < resolution + 0.01 ) { ! m_typeid = final_check( m_corner, atoms, checkfun ); return m_typeid; } ! atomvec_t nbrs; ! ! std::set< int > types; ! atomvec_t::const_iterator ai = atoms.begin(); ! for( ; ai != atoms.end(); ++ai ) { ! int type = checkfun( m_center, m_radius, *ai ); ! ! if( type == TOTAL_INNER ) ! { ! m_typeid = EXCLUDED; ! return EXCLUDED; ! } ! if( type != TOTAL_OUTER ) { ! nbrs.push_back( *ai ); } - - types.insert( type ); - - } ! if( types.count( INNER_SHELL ) == 0 ) ! { ! if( types.count( TOTAL_SHELL ) > 0 ) { ! m_typeid = INCLUDED; return m_typeid; } ! if( types.count( SHELL_OUTER ) == 0 ) { ! m_typeid = EXCLUDED; ! return m_typeid; } } int nincluded = 0; int nexcluded = 0; ! for( int i=0; i < 8; ++i ) { ! int ix = i / 4; ! int iy = ( i % 4 ) / 2; ! int iz = i % 2; ! ! numvec corner(3); ! corner[0] = m_corner[0] + ix * m_length / 2; ! corner[1] = m_corner[1] + iy * m_length / 2; ! corner[2] = m_corner[2] + iz * m_length / 2; ! m_subs[i] = shared_ptr< octree_t >( new octree_t( corner, m_length / 2 ) ); ! ! int subtype = m_subs[i]->divide( nbrs, resolution, checkfun ); ! ! if( subtype == INCLUDED ) { nincluded++; } ! if( subtype == EXCLUDED ) { nexcluded++; } } ! if( nexcluded == 8 ) { m_typeid = EXCLUDED; ! return EXCLUDED; } ! if( nincluded == 8 ) { m_typeid = INCLUDED; ! return INCLUDED; } m_typeid = PARTIAL; return m_typeid; } void octree_t::calculate( potfun_t potfun, std::pair< double, numvec >& min, std::pair< double, numvec >& max ) { ! assert( m_typeid != EXCLUDED ); ! bool nosub = true; ! for( int i=0; i < 8; i++ ) ! { ! if( m_subs[i] != NULL && m_subs[i]->m_typeid != EXCLUDED ) ! { ! nosub = false; ! m_subs[i]->calculate( potfun, min, max ); } } ! if( nosub ) { ! m_potene = potfun( m_corner ); ! if( min.first > m_potene ) ! { ! min.first = m_potene; ! min.second = m_corner; ! } ! if( max.first < m_potene ) ! { ! max.first = m_potene; ! max.second = m_corner; ! } } } ! } // namespace capbox ! } // namespace mort --- 7,437 ---- namespace capbox { + int intpow( int sub, int p ) + { + int r = 1; + for( int i=0; i < p; ++i ) + r *= sub; + return r; + } ! octree_t::octree_t( const numvec& corner, int ndepth, double resolution ) ! : m_corner(corner), m_center(corner), ! m_ndepth(ndepth), m_resolution(resolution) { ! m_ngrid = intpow( 2, m_ndepth ); ! m_length = m_ngrid * resolution; ! m_radius = m_length * sqrt( 0.75 ); ! m_center = m_corner + numvec(m_length/2.0, m_length/2.0, m_length/2.0); ! m_typeid = UNKNOWN_TYPE; } octree_t::~octree_t() { } ! int octree_t::make_shape( shaper_i& shapefun ) { ! // possible change for 2nd cut: ! // excluded -> excluded, ! // included -> included or partial ! // partial -> partial or excluded ! if( m_typeid==EXCLUDED ) { ! // must be excluded ! return EXCLUDED; } ! ! // if m_ndepth is 0, the node should be treated as a point. ! // the result of shapefun can only be INCLUDED or EXCLUDED. ! // PARTIAL is impossible. ! if( m_ndepth==0 ) { ! assert( m_typeid == UNKNOWN_TYPE || m_typeid==INCLUDED ); ! m_typeid = shapefun( m_corner, 0.0 ); ! assert( m_typeid==INCLUDED || m_typeid==EXCLUDED ); return m_typeid; } ! // if m_ndepth > 0, the node is treat as a sphere (with a radius), ! // thus the result of shapefun could be INCLUDED, EXCLUDED or partial ! int type = UNKNOWN_TYPE; ! if( m_typeid==UNKNOWN_TYPE || m_typeid==INCLUDED ) { ! type = shapefun( m_center, m_radius ); ! if( type==EXCLUDED ) { ! // not in the shape, does not matter if 1st or 2nd, ruled out, do not go furthur ! m_typeid = EXCLUDED; ! return m_typeid; } ! if( type==INCLUDED ) { ! m_typeid = type; return m_typeid; } ! assert( type==PARTIAL ); ! ! if( m_typeid==UNKNOWN_TYPE ) { ! m_typeid = PARTIAL; } } + assert( m_typeid==PARTIAL || m_typeid==INCLUDED ); + + + // whenever partial exist, we need to go to the next level. + if( m_subs[0]==NULL ) split(); + int nincluded = 0; int nexcluded = 0; ! for( int i=0; i < 8; ++i ) { ! int subtype = m_subs[i]->make_shape( shapefun ); ! if( subtype==INCLUDED ) { nincluded++; } ! if( subtype==EXCLUDED ) { nexcluded++; } } ! // each node is a cubic box, which is smaller then the sphere, thus ! // it is possible that the box is not really partial, change the result ! // accordingly. ! if( nexcluded==8 ) { m_typeid = EXCLUDED; ! return m_typeid; } ! if( nincluded==8 ) { m_typeid = INCLUDED; ! return m_typeid; } + // now we are sure the node must be partial + // + assert( m_typeid == PARTIAL || m_typeid==INCLUDED ); m_typeid = PARTIAL; return m_typeid; } + + void octree_t::split() + { + assert( m_typeid == INCLUDED || m_typeid==PARTIAL ); + assert( m_ndepth >= 1 && m_ngrid >= 2 ); + for( int i=0; i < 8; ++i ) + { + int ix = i / 4; + int iy = ( i % 4 ) / 2; + int iz = i % 2; + + numvec corner(3); + corner[0] = m_corner[0] + ix * m_length / 2; + corner[1] = m_corner[1] + iy * m_length / 2; + corner[2] = m_corner[2] + iz * m_length / 2; + assert( m_subs[i] == NULL ); + m_subs[i] = octree_ptr( new octree_t( corner, m_ndepth-1, m_resolution ) ); + + // if parent is included (means it is the second shaper), children are included, + // otherwise stay as unknown_type (the init value) + if( m_typeid==INCLUDED ) + { + m_subs[i]->m_typeid = INCLUDED; + } + + if( m_subpol.size() == 0 ) + { + continue; + } + + + // split the value of electro potential + int xstart = ix*m_ngrid/2; + int ystart = iy*m_ngrid/2; + int zstart = iz*m_ngrid/2; + + int subgrid = m_ngrid/2; + m_subs[i]->m_subpol.resize( subgrid*subgrid*subgrid ); + for( int subx=0; subxm_subpol[icrd1] = m_subpol[icrd2]; + } + } + } + } + } + void octree_t::calculate( potfun_t potfun, std::pair< double, numvec >& min, std::pair< double, numvec >& max ) { ! if( m_typeid==EXCLUDED ) ! { ! return; ! } ! ! if( m_typeid==PARTIAL ) ! { ! for( int i=0; i < 8; i++ ) ! { ! if( m_subs[i] != NULL ) ! { ! m_subs[i]->calculate( potfun, min, max ); ! } ! } ! ! return; ! } ! assert( m_typeid==INCLUDED ); ! if( m_subpol.size()==0 ) m_subpol.resize( m_ngrid*m_ngrid*m_ngrid, 0.0 ); ! ! double* psubpol = &m_subpol[0]; ! for( int i=0; i < m_ngrid; ++i) ! { ! double x = m_corner[0] + i*m_resolution; ! for( int j=0; j < m_ngrid; ++j ) ! { ! double y = m_corner[1] + j*m_resolution; ! for( int k=0; k < m_ngrid; ++k ) ! { ! double z = m_corner[2] + k*m_resolution; ! *psubpol = *psubpol + potfun( x, y, z ); ! if( *psubpol < min.first ) ! { ! min.first = *psubpol; ! min.second[0] = x; ! min.second[1] = y; ! min.second[2] = z; ! } ! ! if( max.first < *psubpol ) ! { ! max.first = *psubpol; ! max.second[0] = x; ! max.second[1] = y; ! max.second[2] = z; ! } ! ! psubpol++; ! } } } + + + } + + enum result_e { INNER, PART_INNER, SHELL, PART_SHELL, OUTER }; + + shell_shaper::shell_shaper( const double* pcrd, const double* prad, int nsphere, double ion_radius, double shell_extent ) + : m_exclusion( nsphere, 0 ) + { + m_psphere = pcrd; + m_pradius = prad; + + m_nsphere = nsphere; + + m_ion_radius = ion_radius; + + double max_radius = *std::max_element( prad, prad+nsphere ); + m_outer_radius = max_radius + ion_radius + shell_extent; + } + + int shell_shaper::operator()( const numvec& center, double radius ) + { + vector num(5, 0); + + // the exclusion list cotains the ID of spheres who are totally + // out of the previous node. If the previous node happen to be + // parent of the current node, its exclusion list should be + // inherited since since children is always inside the parent. + + // go back until find my father, whose radius is bigger than mine + // or there might not be one. + int old_depth = m_exclusion_radii.size(); + while( m_exclusion_radii.size()>0 && radius>=m_exclusion_radii.back() ) + { + m_exclusion_radii.pop_back(); + } + + + // clear the exclusion of my brothers and their descents. leave only + // the exclusion of my father and his parents + int new_depth = m_exclusion_radii.size(); + if( new_depth != old_depth ) + { + assert( new_depth < old_depth ); + assert( m_exclusion.size()==m_nsphere ); + + for(int is=0; is < m_nsphere; ++is ) + { + if( m_exclusion[is] > new_depth ) + m_exclusion[is] = 0; + } + } + + + vector own_exclusion; + for( int i=0; i < m_nsphere; ++i ) + { + // if i is in my father's exclusion list, skip it + if( m_exclusion[i]> 0 ) + { + continue; + } + + + int type = check_each( center, radius, i ); ! // if I am sucked inside any sphere, I must be excluded ! if( type == INNER ) { ! return EXCLUDED; ! } ! ! num[type] += 1; ! ! ! // if I am a point, no need to track exclusion list ! // because I have no children. ! if( radius >0.0 && type == OUTER ) ! own_exclusion.push_back( i ); ! } ! ! if( own_exclusion.size()>0 ) ! { ! m_exclusion_radii.push_back( radius ); ! ! ! int ndepth = m_exclusion_radii.size(); ! ! for( int i=0; i < own_exclusion.size(); ++i ) ! { ! int is = own_exclusion[i]; ! // if a sphere is in my own exclusion list, it cannot be in my father's exclusion ! // list, since I will skip it. ! assert( m_exclusion[is]==0 ); ! m_exclusion[is] = ndepth; } + + /* + for( int i=0; i < m_exclusion_radii.size(); ++i ) std::cout << " "; + std::cout << "x,y,z,r, self_exclusion, total_exclusion:"; + std::cout << center[0] << " " << center[1] << " " << center[2] << " " << radius << " "; + std::cout << my_exclusion.size() << " " << m_exclusion.size() << std::endl; */ } + + + + assert( num[INNER]==0 ); + + if( radius == 0.0 ) + { + assert( num[PART_INNER]==0 && num[PART_SHELL]==0 ); + // only two possiblities is SHELL and OUTER, + // if nshell==0, all outer, then exclude + + return num[SHELL]==0 ? EXCLUDED : INCLUDED; + } + + if( num[PART_INNER] > 0 ) + { + // intersect with an atom, must be partial. + return PARTIAL; + } + + if( num[SHELL] > 0 ) + { + // inside some atoms' shell, must be included? taken from old leap code + return INCLUDED; + } ! // only possible: PART_SHELL, OUTER, if no part_shell, all outer ! // should be excluded ! return (num[PART_SHELL]==0) ? EXCLUDED : PARTIAL; ! } ! ! int shell_shaper::check_each( const numvec& center, double radius, int isphere ) const ! { ! double dx = center[0] - m_psphere[3*isphere]; ! double dy = center[1] - m_psphere[3*isphere+1]; ! double dz = center[2] - m_psphere[3*isphere+2]; ! double d = sqrt( dx*dx + dy*dy + dz*dz ); ! ! double inner_radius = m_pradius[isphere] + m_ion_radius; ! ! ! if( radius==0.0 ) ! { ! if( d < inner_radius ) ! { ! return INNER; ! } ! ! if( d < m_outer_radius ) return SHELL; ! ! return OUTER; ! } ! ! if( d + radius < inner_radius ) return INNER; ! ! if( d - radius < inner_radius ) return PART_INNER; ! ! if( d + radius < m_outer_radius ) return SHELL; ! ! if( d - radius < m_outer_radius ) return PART_SHELL; ! ! return OUTER; ! } + outsphere_shaper::outsphere_shaper( const numvec& sphere, double radius ) + :m_sphere( sphere ), m_radius(radius) + { + } + int outsphere_shaper::operator()( const numvec& center, double radius ) + { + double dx = center[0] - m_sphere[0]; + double dy = center[1] - m_sphere[1]; + double dz = center[2] - m_sphere[2]; + double d = sqrt( dx*dx + dy*dy + dz*dz ); + + if( radius==0.0 ) + { + return (d potfun_t; class octree_t { public: ! octree_t( const numvec& corner, double length ); ! ~octree_t(); ! int divide( const atomvec_t& atoms, double resolution, const checkfun_t& check ); ! void calculate( potfun_t potfun, std::pair< double, numvec >& min, std::pair< double, numvec >& max ); ! std::pair< double, numvec > minpot(); ! ! std::pair< double, numvec > maxpot(); private: --- 8,30 ---- namespace capbox { ! class shaper_i; ! typedef function< double ( double, double, double ) > potfun_t; class octree_t { public: ! octree_t( const numvec& corner, int ndepth, double resolution ); ! virtual ~octree_t(); ! void split( ); ! int make_shape( shaper_i& shapefun ); ! void calculate( potfun_t potfun, std::pair< double, numvec >& min, std::pair< double, numvec >& max ); private: *************** namespace mort *** 105,147 **** numvec m_corner; numvec m_center; ! ! atomvec_t m_atoms; ! double m_length; double m_radius; ! double m_potene; - int m_typeid; }; struct elepot_t { ! elepot_t( const molecule_t& mol ) { ! copy( mol.atom_begin(), mol.atom_end(), back_inserter( m_atoms ) ); } ! elepot_t( const moref_t& resd ) { ! assert( resd.cmpid() == RESD ); ! copy( resd.related_atom_begin(), resd.related_atom_end(), back_inserter( m_atoms ) ); } ! double operator()( const numvec& pos ) { double elep = 0.0; ! for( int i=0; i < (int)m_atoms.size(); ++i ) { ! double dist = mort::dist( pos, m_atoms[i].get_v(POSITION) ); ! elep = elep + m_atoms[i].get_d(PCHG) / (dist*dist); } return elep; } ! atomvec_t m_atoms; }; } // namespace capbox --- 33,109 ---- numvec m_corner; numvec m_center; ! ! int m_typeid; ! ! int m_ndepth; ! ! int m_ngrid; ! double m_length; double m_radius; + + double m_resolution; ! vector m_subpol; }; + + typedef shared_ptr< octree_t > octree_ptr; struct elepot_t { ! elepot_t( const molecule_t& mol, int natom=0 ) ! : m_natom(natom) { ! if( m_natom==0 ) ! { ! m_natom = mol.natom(); ! } ! ! m_pchg = get_dptr( mol, ATOM, PCHG ); ! m_pcrd = get_vptr( mol, ATOM, POSITION ); } ! elepot_t( const numvec& pos, double chg ) ! : m_spos( pos ), m_schg(chg) { ! m_natom = 0; } ! double operator()( double x, double y, double z ) { + if( m_natom==0 ) + { + double dx = m_spos[0] - x; + double dy = m_spos[1] - y; + double dz = m_spos[2] - z; + double d2 = dx*dx + dy*dy + dz*dz; + return m_schg/d2; + } + + double elep = 0.0; ! for( int i=0; i < m_natom; ++i ) { ! double dx = m_pcrd[3*i ] - x; ! double dy = m_pcrd[3*i+1] - y; ! double dz = m_pcrd[3*i+2] - z; ! double d2 = dx*dx + dy*dy + dz*dz; ! elep = elep + m_pchg[i] / d2; } + return elep; } ! int m_natom; ! const double* m_pchg; ! const double* m_pcrd; ! ! // this is for single atom: 's' stands for single. ! numvec m_spos; ! double m_schg; }; } // namespace capbox *************** namespace mort *** 150,152 **** --- 112,185 ---- #endif + #include + + namespace mort + { + namespace capbox + { + + enum shaper_result_e { UNKNOWN_TYPE, INCLUDED, PARTIAL, EXCLUDED }; + + class shaper_i + { + public: + + shaper_i() {} + + virtual ~shaper_i() {} + + virtual int operator()( const numvec& center, double radius ) = 0; + }; + + + class shell_shaper : public shaper_i + { + public: + + shell_shaper( const double* pcrd, const double* prad, int nsphere, double ion_radius, double shell_extent ); + + virtual ~shell_shaper() {} + + virtual int operator()( const numvec& center, double radius ); + + private: + + int check_each( const numvec& center, double radius, int isphere ) const; + + const double* m_psphere; + + const double* m_pradius; + + int m_nsphere; + + double m_ion_radius; + + double m_outer_radius; + + vector m_exclusion; + + vector m_exclusion_radii; + + }; + + class outsphere_shaper : public shaper_i + { + public: + + outsphere_shaper( const numvec& sphere, double radius ); + + virtual ~outsphere_shaper() {} + + virtual int operator()( const numvec& center, double radius ); + + private: + + numvec m_sphere; + + double m_radius; + }; + + } // namespace capbox + + } // namespace mort Index: mortsrc/capbox/solvate.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/capbox/solvate.cpp,v retrieving revision 1.17 diff -p -r1.17 solvate.cpp *** src/gleap/mortsrc/capbox/solvate.cpp 9 Apr 2008 23:34:00 -0000 1.17 --- src/gleap/mortsrc/capbox/solvate.cpp 29 Apr 2008 19:20:55 -0000 *************** *** 6,77 **** namespace mort { - using std::logic_error; - using boost::format; ! double get_parm99_vdwr(const moref_t& a) { ! static const char types[][4] = { "H" , "HO", "HS", "HC", ! "H1", "H2", "H3", "HP", ! "HA", "H4", "H5", "HW", ! "HZ", "O", "O2", "OW", ! "OH", "OS" }; ! ! ! static const double radiis[] = { 0.6, 1.0, 0.6, 1.487, ! 1.387, 1.287, 1.187, 1.10, ! 1.459, 1.409, 1.359, 1.0, ! 1.459, 1.6612,1.6612,1.7683, ! 1.721, 1.6837}; ! ! static const int size = sizeof( radiis ) / sizeof(double); ! ! assert( size == 18 ); ! ! string type = a.get_s(TYPE); ! ! if( type[0]=='C' || type[0]=='c' ) return 1.9080; ! ! if( type[0]=='N' || type[0]=='n' ) return 1.8240; ! ! if( type[0]=='S' || type[0]=='s' ) return 2.0; ! ! if( type[0]=='P' || type[0]=='p' ) return 2.1; ! ! if( type[0]=='o' ) return 1.7; ! ! if( type[0]=='h' ) return 1.0; ! ! int pos = find( types, types+size, type ) - types; ! ! if( pos == size ) ! { ! int elem = a.get_i(ELEMENT); ! if( elem==CARBON ) return 1.9080; ! ! if( elem==NITROGEN ) return 1.8240; ! ! if( elem==SULFUR ) return 2.0; ! ! if( elem==15 ) return 2.1; ! ! if( elem==OXYGEN ) return 1.7; ! if( elem==HYDROGEN ) return 1.0; ! throw logic_error( "Error: unknown vdwr for type " + type ); ! } ! return radiis[pos]; } ! void set_parm99_vdwr( molecule_t& m ) { ! moiter_t a = m.atom_begin(); ! for( ; a != m.atom_end(); ++a ) { ! a->set_d(VDWR, get_parm99_vdwr(*a) ); } } --- 6,51 ---- namespace mort { using boost::format; ! double get_amberff_vdwr(const string& type) { ! molecule_ptr pff = mortenv().get_mol("_amber-atom"); ! if( pff==NULL ) ! { ! std::cout << "Warning: no force field has been loaded, VDW radii set to 1.5." << std::endl; ! return 1.5; ! } ! moref_t pa( *pff, ATOM, -1 ); ! if( pff->get_atom( type, pa ) ) ! { ! // in parm99.dat, some atom types has zero radius (HO and HW for example), which would be too small ! // change them to 1.0 here. ! double r = pa.get_d(RSTAR); ! return ( r<0.01) ? 1.0 : r; ! } ! // for extra points ! if( type[0]=='E' ) ! { ! return 0.2; ! } ! std::cout << "Warning: unknown amber atom type " << type << ", using radii 1.5" << std::endl; ! return 1.5; } ! void set_amberff_vdwr( molecule_t& m, int start ) { ! m.atom_begin()->set_d(VDWR, 0.0); ! ! double* pvdwr = get_dptr( m, ATOM, VDWR ); ! const string* ptype = get_sptr( m, ATOM, TYPE ); ! int natom = m.natom(); ! for( int i=start; i < natom; ++i ) { ! pvdwr[i] = get_amberff_vdwr(ptype[i]); } } *************** namespace mort *** 92,103 **** corner[2] = center[2] + 0.5*box[2]*nz; check_shape shaper( shape, center, extent ); ! std::cout << "nx,ny,nz:"; ! std::cout << format( "%8d" ) % nx; ! std::cout << format( "%8d" ) % ny; ! std::cout << format( "%8d" ) % nz; ! std::cout << std::endl; ! int natom = mol.natom(); for( int ix=0; ix <= nx; ix++ ) --- 66,72 ---- corner[2] = center[2] + 0.5*box[2]*nz; check_shape shaper( shape, center, extent ); ! // std::cout << "nx,ny,nz:" << format( "%8d %8d %8d" ) % nx % ny % nz << std::endl; int natom = mol.natom(); for( int ix=0; ix <= nx; ix++ ) *************** namespace mort *** 111,123 **** cnt[1] = corner[1] - (iy+0.5)*box[1]; cnt[2] = corner[2] - (iz+0.5)*box[2]; ! /* ! std::cout << "place solvent on center: "; ! std::cout << format( "%8.3f" ) % cnt[0]; ! std::cout << format( "%8.3f" ) % cnt[1]; ! std::cout << format( "%8.3f" ) % cnt[2]; ! std::cout << std::endl; ! */ check_collision collide( mol, natom, shell_extent, closeness, cnt, box ); add_solvent( mol, solvent, cnt, shaper, collide ); --- 80,87 ---- cnt[1] = corner[1] - (iy+0.5)*box[1]; cnt[2] = corner[2] - (iz+0.5)*box[2]; ! ! //std::cout << "place solvent on center: %8.3f %8.3f %8.3f" % cnt[0] % cnt[1] % cnt[2] << std::endl; check_collision collide( mol, natom, shell_extent, closeness, cnt, box ); add_solvent( mol, solvent, cnt, shaper, collide ); *************** namespace mort *** 128,160 **** void add_solvent( molecule_t& mol, const molecule_t& solvent, const numvec& cnt, check_shape shaper, check_collision collide ) { ! numvec dis = cnt - center( solvent ); const double* solventpos = get_vptr(solvent, ATOM, POSITION); ! moiter_t resd = solvent.resd_begin(); ! for( ; resd != solvent.resd_end(); ++resd ) { bool add = true; ! moiter_t atom = resd->related_atom_begin(); ! ! for( ; atom != resd->related_atom_end(); ++atom ) { ! int pid = atom->absid(); int pid3 = pid*3; ! double curt[3]; ; curt[0] = dis[0]+solventpos[pid3]; curt[1] = dis[1]+solventpos[pid3+1]; curt[2] = dis[2]+solventpos[pid3+2]; ! if( ! shaper( curt, get_parm99_vdwr(*atom) ) ) { add = false; break; } ! if( ! collide( curt, get_parm99_vdwr(*atom) ) ) { add = false; break; --- 92,124 ---- void add_solvent( molecule_t& mol, const molecule_t& solvent, const numvec& cnt, check_shape shaper, check_collision collide ) { ! numvec dis = cnt - center(solvent); const double* solventpos = get_vptr(solvent, ATOM, POSITION); ! moiter_t ri = solvent.resd_begin(); ! for( ; ri != solvent.resd_end(); ++ri ) { bool add = true; ! moiter_t ai = ri->related_atom_begin(); ! for( ; ai != ri->related_atom_end(); ++ai ) { ! int pid = ai->absid(); int pid3 = pid*3; ! double curt[3]; curt[0] = dis[0]+solventpos[pid3]; curt[1] = dis[1]+solventpos[pid3+1]; curt[2] = dis[2]+solventpos[pid3+2]; ! double solvent_vdwr = get_amberff_vdwr( ai->get_s(TYPE) ); ! if( ! shaper(curt, solvent_vdwr) ) { add = false; break; } ! if( ! collide(curt, solvent_vdwr) ) { add = false; break; *************** namespace mort *** 163,169 **** if( add ) { ! moref_t copy = merge( mol, *resd ); move( copy, dis ); } } --- 127,133 ---- if( add ) { ! moref_t copy = merge( mol, *ri ); move( copy, dis ); } } *************** namespace mort *** 240,252 **** m_closeness( closeness ) { double max_solute_r = 0.0; ! ! moiter_t a = mol.atom_begin(); ! for( int i=0; i < natom; ++i, ++a ) { ! if( a->get_d(VDWR)*m_closeness > max_solute_r ) { ! max_solute_r = a->get_d(VDWR)*m_closeness; } } --- 204,215 ---- m_closeness( closeness ) { double max_solute_r = 0.0; ! const double* pr = get_dptr( mol, ATOM, VDWR ); ! for( int i=0; i < natom; ++i ) { ! if( pr[i]*m_closeness > max_solute_r ) { ! max_solute_r = pr[i]*m_closeness; } } *************** namespace mort *** 264,272 **** double zmax = center[2] + 0.5*box[2] + buffer; const double* pos = get_vptr( mol, ATOM, POSITION ); ! ! a = mol.atom_begin(); ! for( int i=0; i < natom; ++i,++a ) { int i3 = 3*i; if( pos[i3 ]xmax ) continue; --- 227,234 ---- double zmax = center[2] + 0.5*box[2] + buffer; const double* pos = get_vptr( mol, ATOM, POSITION ); ! ! for( int i=0; i < natom; ++i ) { int i3 = 3*i; if( pos[i3 ]xmax ) continue; *************** namespace mort *** 276,282 **** m_solute_pos.push_back( pos[i3 ] ); m_solute_pos.push_back( pos[i3+1] ); m_solute_pos.push_back( pos[i3+2] ); ! m_solute_vdwr.push_back( a->get_d(VDWR)*m_closeness ); } } --- 238,244 ---- m_solute_pos.push_back( pos[i3 ] ); m_solute_pos.push_back( pos[i3+1] ); m_solute_pos.push_back( pos[i3+2] ); ! m_solute_vdwr.push_back( pr[i]*m_closeness ); } } *************** namespace mort *** 410,444 **** void solvate_cap( molecule_t& mol, const molecule_t& solvent, const numvec& center, double radius, double closeness ) { numvec size( 2*radius, 2*radius, 2*radius ); capbox::solvate( mol, solvent, center, size, CAP, 0.0, closeness ); mol.set_i(SOLUTE, CAP); mol.set_v(CAP, numvec(center[0], center[1], center[2], radius) ); } void solvate_box( molecule_t& mol, const molecule_t& solvent, double buffer, double closeness ) { ! set_parm99_vdwr( mol ); ! numvec extent = center_by_vdwr( mol ); ! extent[0] += 2*buffer; ! extent[1] += 2*buffer; ! extent[2] += 2*buffer; capbox::solvate( mol, solvent, numvec(0.0,0.0,0.0), extent, BOX, 0.0, closeness ); mol.set_i(SOLUTE, BOX); mol.set_v(BOX, numvec( extent[0], extent[1], extent[2], 90.0 ) ); ! set_parm99_vdwr( mol ); center_by_vdwr( mol ); } void solvate_oct( molecule_t& mol, const molecule_t& solvent, double buffer, double closeness ) { ! set_parm99_vdwr( mol ); numvec cnt = center( mol ); move( mol, -1.0 * cnt ); --- 372,418 ---- void solvate_cap( molecule_t& mol, const molecule_t& solvent, const numvec& center, double radius, double closeness ) { + int solute_size = mol.natom(); + mol.set_i( SOLUTESIZE, solute_size ); + + set_amberff_vdwr( mol ); + numvec size( 2*radius, 2*radius, 2*radius ); capbox::solvate( mol, solvent, center, size, CAP, 0.0, closeness ); + // set vdw radius again for newly added solvent molecules. + set_amberff_vdwr( mol, solute_size ); + mol.set_i(SOLUTE, CAP); mol.set_v(CAP, numvec(center[0], center[1], center[2], radius) ); } void solvate_box( molecule_t& mol, const molecule_t& solvent, double buffer, double closeness ) { ! int solute_size = mol.natom(); ! mol.set_i( SOLUTESIZE, solute_size ); ! set_amberff_vdwr( mol ); ! ! numvec extent = center_by_vdwr(mol); ! extent += numvec( 2*buffer, 2*buffer, 2*buffer ); capbox::solvate( mol, solvent, numvec(0.0,0.0,0.0), extent, BOX, 0.0, closeness ); mol.set_i(SOLUTE, BOX); mol.set_v(BOX, numvec( extent[0], extent[1], extent[2], 90.0 ) ); ! set_amberff_vdwr( mol ); center_by_vdwr( mol ); } void solvate_oct( molecule_t& mol, const molecule_t& solvent, double buffer, double closeness ) { ! int solute_size = mol.natom(); ! mol.set_i( SOLUTESIZE, solute_size ); ! ! set_amberff_vdwr( mol ); numvec cnt = center( mol ); move( mol, -1.0 * cnt ); *************** namespace mort *** 479,503 **** capbox::solvate( mol, solvent, numvec(0.0, 0.0, 0.0), extent, OCT, 0.0, closeness ); ! set_parm99_vdwr( mol ); center_by_vdwr( mol ); ewald_rotate( mol ); mol.set_i(SOLUTE, BOX); double boxlen = extmax * sqrt(0.75); mol.set_v(BOX, numvec(boxlen, boxlen, boxlen, 109.471219) ); } void solvate_shell( molecule_t& mol, const molecule_t& solvent, double buffer, double closeness ) { ! set_parm99_vdwr( mol ); numvec extent = center_by_vdwr( mol ); extent += numvec( 2*buffer, 2*buffer, 2*buffer ); capbox::solvate( mol, solvent, numvec(0.0, 0.0, 0.0), extent, BOX, buffer, closeness ); ! set_parm99_vdwr( mol ); center_by_vdwr( mol ); } --- 453,483 ---- capbox::solvate( mol, solvent, numvec(0.0, 0.0, 0.0), extent, OCT, 0.0, closeness ); ! // assign vdwr for solvent atoms and center again. ! set_amberff_vdwr( mol, solute_size ); center_by_vdwr( mol ); ewald_rotate( mol ); mol.set_i(SOLUTE, BOX); + + double boxlen = extmax * sqrt(0.75); mol.set_v(BOX, numvec(boxlen, boxlen, boxlen, 109.471219) ); } void solvate_shell( molecule_t& mol, const molecule_t& solvent, double buffer, double closeness ) { ! int solute_size = mol.natom(); ! mol.set_i( SOLUTESIZE, solute_size ); ! ! set_amberff_vdwr( mol ); numvec extent = center_by_vdwr( mol ); extent += numvec( 2*buffer, 2*buffer, 2*buffer ); capbox::solvate( mol, solvent, numvec(0.0, 0.0, 0.0), extent, BOX, buffer, closeness ); ! set_amberff_vdwr( mol ); center_by_vdwr( mol ); } Index: mortsrc/capbox/solvate.hpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/capbox/solvate.hpp,v retrieving revision 1.5 diff -p -r1.5 solvate.hpp *** src/gleap/mortsrc/capbox/solvate.hpp 9 Apr 2008 23:34:00 -0000 1.5 --- src/gleap/mortsrc/capbox/solvate.hpp 29 Apr 2008 19:20:55 -0000 *************** namespace mort *** 60,71 **** --- 60,76 ---- double m_closeness; }; + void solvate( molecule_t& mol, const molecule_t& solvent, const numvec& center, const numvec& length, int shape, double shell_extent, double closeness ); void add_solvent( molecule_t& mol, const molecule_t& solvent, const numvec& pos, check_shape shaper, check_collision collide ); } // namespace capbox + void set_amberff_vdwr( molecule_t& m, int start=0 ); + + double get_amberff_vdwr( const string& type ); + } // namespace mort #endif Index: mortsrc/common/hash.txt =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/common/hash.txt,v retrieving revision 1.6 diff -p -r1.6 hash.txt *** src/gleap/mortsrc/common/hash.txt 21 Jan 2008 22:20:10 -0000 1.6 --- src/gleap/mortsrc/common/hash.txt 29 Apr 2008 19:20:55 -0000 *************** shell *** 96,101 **** --- 96,102 ---- smarts smiles solute + solutesize stereo stick strbnd Index: mortsrc/common/hashcode.hpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/common/hashcode.hpp,v retrieving revision 1.10 diff -p -r1.10 hashcode.hpp *** src/gleap/mortsrc/common/hashcode.hpp 21 Jan 2008 22:20:10 -0000 1.10 --- src/gleap/mortsrc/common/hashcode.hpp 29 Apr 2008 19:20:55 -0000 *************** namespace mort *** 107,112 **** --- 107,113 ---- static const hashid_t SMARTS = 1121097474LLU; static const hashid_t SMILES = 1095633666LLU; static const hashid_t SOLUTE = 274725306LLU; + static const hashid_t SOLUTESIZE = 477433944480186LLU; static const hashid_t STEREO = 854043966LLU; static const hashid_t STICK = 16900542LLU; static const hashid_t STRBND = 203302926LLU; Index: mortsrc/guilib/mainwin.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/guilib/mainwin.cpp,v retrieving revision 1.11 diff -p -r1.11 mainwin.cpp *** src/gleap/mortsrc/guilib/mainwin.cpp 2 Apr 2008 18:25:11 -0000 1.11 --- src/gleap/mortsrc/guilib/mainwin.cpp 29 Apr 2008 19:20:55 -0000 *************** namespace mort *** 14,21 **** database_t& content() { ! static database_t g_content; ! return g_content; } string find_file( const string& name ) --- 14,22 ---- database_t& content() { ! //static database_t g_content; ! //return g_content; ! return mortenv(); } string find_file( const string& name ) Index: mortsrc/object/component.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/object/component.cpp,v retrieving revision 1.15 diff -p -r1.15 component.cpp *** src/gleap/mortsrc/object/component.cpp 31 Aug 2007 21:44:50 -0000 1.15 --- src/gleap/mortsrc/object/component.cpp 29 Apr 2008 19:20:55 -0000 *************** namespace mort *** 89,95 **** void compress_byid_v( vector& src, const vector& ids, int length ) { ! assert( length>0 && src.size()>length*ids.size() ); vector< double > dst( length*ids.size() ); --- 89,95 ---- void compress_byid_v( vector& src, const vector& ids, int length ) { ! assert( length>0 && src.size()>=length*ids.size() ); vector< double > dst( length*ids.size() ); *************** namespace mort *** 173,188 **** return m_absids.end()-n; } ! vector::iterator component_t::insert( const vector::iterator& pos, int n ) { ! int seq = pos - m_absids.begin(); ! m_absids.insert( pos, n, m_idcounter ); ! vector::iterator j = m_absids.begin()+seq; ! for(int i=0; i < n; ++i,++j) { *j = m_idcounter; ++m_idcounter; } --- 173,196 ---- return m_absids.end()-n; } ! vector::iterator component_t::insert( int seq, int n ) { ! vector::iterator oldpos = m_absids.begin()+seq; ! m_absids.insert( oldpos, n, -1 ); ! // it is possible that the memory has been reallocated ! // and moved to new address, thus it is necessary to create ! // nwe pos ! ! vector::iterator newpos = m_absids.begin() + seq; ! ! ! // assign new absid ! vector::iterator j = newpos; for(int i=0; i < n; ++i,++j) { + assert( *j == -1 ); *j = m_idcounter; ++m_idcounter; } *************** namespace mort *** 190,196 **** resize_parm(); m_clean = false; ! return j; } bool component_t::remove( int id ) --- 198,204 ---- resize_parm(); m_clean = false; ! return newpos; } bool component_t::remove( int id ) *************** namespace mort *** 199,205 **** if( i == m_absids.end() ) return false; m_absids.erase( i ); - m_clean = false; return true; } --- 207,212 ---- *************** namespace mort *** 233,243 **** void component_t::cleanup(map& old2new) { - if( m_clean ) - { - return; - } - compress_content_n( m_icontent, m_absids ); compress_content_n( m_dcontent, m_absids ); compress_content_n( m_scontent, m_absids ); --- 240,245 ---- *************** namespace mort *** 259,265 **** adjust_occupy( m_voccupy, old2new ); adjust_occupy( m_aoccupy, old2new ); - m_clean = true; } void component_t::set_i(const hashid_t& parmid, int absid, const int& value) --- 261,266 ---- *************** namespace mort *** 545,550 **** --- 546,560 ---- return pv; } + string* component_t::get_sptr(const hashid_t& parmid, int absid) + { + string* ps = NULL; + if( !get_sptr(parmid, absid, ps) ) + throw logic_error( "Error: undefinded numvec parameter " + unhash(parmid) ); + return ps; + } + + double const* component_t::get_dptr(const hashid_t& parmid, int absid) const { const double* pd = NULL; *************** namespace mort *** 553,558 **** --- 563,570 ---- return pd; } + + const double* component_t::get_vptr(const hashid_t& parmid, int absid) const { const double* pv = NULL; *************** namespace mort *** 560,565 **** --- 572,587 ---- throw logic_error( "Error: undefinded numvec parameter " + unhash(parmid) ); return pv; } + + string const* component_t::get_sptr(const hashid_t& parmid, int absid) const + { + const string* ps = NULL; + if( !get_sptr(parmid, absid, ps) ) + throw logic_error( "Error: undefinded numvec parameter " + unhash(parmid) ); + return ps; + } + + const any& component_t::get_a(const hashid_t& parmid, int absid) const { Index: mortsrc/object/component.hpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/object/component.hpp,v retrieving revision 1.8 diff -p -r1.8 component.hpp *** src/gleap/mortsrc/object/component.hpp 4 May 2007 12:18:04 -0000 1.8 --- src/gleap/mortsrc/object/component.hpp 29 Apr 2008 19:20:55 -0000 *************** namespace mort *** 54,64 **** void swap( component_t& rhs ); /// \brief add new object to the end. ! iterator append(int n); /// \brief insert new object at the certain position. /// \param i the position where new object will be added. ! iterator insert(const iterator& pos, int n); /// \brief remove object with given ID number /// \param the ID number of the object going to be removed. --- 54,64 ---- void swap( component_t& rhs ); /// \brief add new object to the end. ! iterator append(int n=1); /// \brief insert new object at the certain position. /// \param i the position where new object will be added. ! iterator insert(int relid, int n=1); /// \brief remove object with given ID number /// \param the ID number of the object going to be removed. *************** namespace mort *** 117,124 **** --- 117,127 ---- double* get_dptr( const hashid_t& parmid, int absid ); double* get_vptr( const hashid_t& parmid, int absid ); + string* get_sptr( const hashid_t& parmid, int absid ); + double const* get_dptr(const hashid_t& parmid, int absid) const; double const* get_vptr(const hashid_t& parmid, int absid) const; + string const* get_sptr(const hashid_t& parmid, int absid) const; // Testing functions bool has_i( const hashid_t& parmid, int absid ) const; Index: mortsrc/object/database.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/object/database.cpp,v retrieving revision 1.12 diff -p -r1.12 database.cpp *** src/gleap/mortsrc/object/database.cpp 10 Feb 2008 15:04:43 -0000 1.12 --- src/gleap/mortsrc/object/database.cpp 29 Apr 2008 19:20:55 -0000 *************** namespace mort *** 153,159 **** } moref_t resd = get_resd( mask.substr( 0, dot ) ); ! string name = mask.substr( dot + 1, mask.length() - dot - 1 ); return resd.get_related_atom( sparm_comper1(NAME, name) ); --- 153,159 ---- } moref_t resd = get_resd( mask.substr( 0, dot ) ); ! string name = mask.substr( dot + 1, mask.length() - dot - 1 ); return resd.get_related_atom( sparm_comper1(NAME, name) ); *************** namespace mort *** 177,185 **** throw logic_error("Error: can not find molecule " + molname + " in database" ); } ! string rid = mask.substr( dot + 1, mask.length() - dot - 1 ); ! return moref_t( *pmol, RESD, atoi( rid.c_str() ) - 1 ); } void database_t::set( const string& name, const entity_ptr& pe ) --- 177,193 ---- throw logic_error("Error: can not find molecule " + molname + " in database" ); } ! string r = mask.substr( dot + 1, mask.length() - dot - 1 ); ! assert( r.length()>0 ); ! ! if( isdigit(r[0]) ) ! { ! int rid = atoi( r.c_str() ); ! assert( rid > 0 ); ! return moref_t( *pmol, RESD, rid-1 ); ! } ! return pmol->get_resd( r ); } void database_t::set( const string& name, const entity_ptr& pe ) *************** namespace mort *** 271,276 **** --- 279,290 ---- return bnds.size()>0; } + + database_t& mortenv() + { + static database_t g_mortenv; + return g_mortenv; + } Index: mortsrc/object/database.hpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/object/database.hpp,v retrieving revision 1.14 diff -p -r1.14 database.hpp *** src/gleap/mortsrc/object/database.hpp 17 Aug 2007 03:04:43 -0000 1.14 --- src/gleap/mortsrc/object/database.hpp 29 Apr 2008 19:20:55 -0000 *************** namespace mort *** 109,114 **** --- 109,115 ---- typedef shared_ptr database_ptr; + database_t& mortenv(); void transform( database_t& mdb, const double* mat ); Index: mortsrc/object/entity.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/object/entity.cpp,v retrieving revision 1.9 diff -p -r1.9 entity.cpp *** src/gleap/mortsrc/object/entity.cpp 30 Aug 2007 03:23:02 -0000 1.9 --- src/gleap/mortsrc/object/entity.cpp 29 Apr 2008 19:20:55 -0000 *************** namespace mort *** 547,558 **** return m_acontent.end(); } - entity_t& mortenv() - { - static entity_t g_mortenv; - return g_mortenv; - } - bool ismol( const entity_ptr& pe ) { molecule_ptr pm = dynamic_pointer_cast< molecule_t >( pe ); --- 547,552 ---- Index: mortsrc/object/entity.hpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/object/entity.hpp,v retrieving revision 1.6 diff -p -r1.6 entity.hpp *** src/gleap/mortsrc/object/entity.hpp 30 Aug 2007 03:22:49 -0000 1.6 --- src/gleap/mortsrc/object/entity.hpp 29 Apr 2008 19:20:55 -0000 *************** *** 4,9 **** --- 4,10 ---- #include #include #include + #include #include #include #include *************** namespace mort *** 14,19 **** --- 15,21 ---- { using std::map; using std::string; + using std::ostream; using boost::any; using boost::shared_ptr; *************** namespace mort *** 172,181 **** map m_acontent; }; ! entity_t& mortenv(); ! ! typedef shared_ptr entity_ptr; bool ismol( const entity_ptr& pe ); --- 174,182 ---- map m_acontent; }; ! ostream& mortlog( int level ); + typedef shared_ptr entity_ptr; bool ismol( const entity_ptr& pe ); Index: mortsrc/object/molecule.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/object/molecule.cpp,v retrieving revision 1.15 diff -p -r1.15 molecule.cpp *** src/gleap/mortsrc/object/molecule.cpp 20 Aug 2007 18:22:40 -0000 1.15 --- src/gleap/mortsrc/object/molecule.cpp 29 Apr 2008 19:20:55 -0000 *************** namespace mort *** 349,359 **** return *moiter_t(*this, PTOR, getcmp(PTOR)->append(1) ); } void molecule_t::remove(const moref_t& mo) { int cmpid = mo.cmpid(); int absid = mo.absid(); ! getcmp(cmpid)->remove(absid); const set* pnbr = NULL; --- 349,383 ---- return *moiter_t(*this, PTOR, getcmp(PTOR)->append(1) ); } + moref_t molecule_t::insert_atom( int relid ) + { + assert( relid>0 && relid < natom() ); + + vector::iterator i = getcmp(ATOM)->insert(relid); + assert( i - getcmp(ATOM)->absid_begin()==relid ); + // assert( *i == natom()-1 ); + + return moref_t(*this, ATOM, *i); + } + + + moref_t molecule_t::insert_resd( int relid ) + { + assert( relid>0 && relid < nresd() ); + + vector::iterator i = getcmp(RESD)->insert(relid); + + assert( i - getcmp(RESD)->absid_begin()==relid ); + // assert( *i == nresd()-1 ); + + return moref_t(*this, RESD, *i); + } + void molecule_t::remove(const moref_t& mo) { int cmpid = mo.cmpid(); int absid = mo.absid(); ! getcmp(cmpid)->remove(absid); const set* pnbr = NULL; *************** namespace mort *** 397,410 **** remove(a); } void molecule_t::remove_resd(const moref_t& r) { assert( r.cmpid() == RESD ); vector ratoms( r.related_atom_begin(), r.related_atom_end() ); ! std::for_each( ratoms.begin(), ratoms.end(), bind( mem_fn(&molecule_t::remove), this, _1 ) ); ! remove(r); } --- 421,453 ---- remove(a); } + + void molecule_t::remove_bond(const moref_t& b) + { + adjacency_t* b2a = getadj( BOND, ATOM ); + vector* pas = b2a->getvec( b.absid() ); + assert( pas->size()==2 ); + int atom_1 = pas->at( 0 ); + int atom_2 = pas->at( 1 ); + + + adjacency_t* a2a = getadj( ATOM, ATOM ); + a2a->remove( atom_1, atom_2 ); + a2a->remove( atom_2, atom_1 ); + + remove( b ); + } void molecule_t::remove_resd(const moref_t& r) { assert( r.cmpid() == RESD ); vector ratoms( r.related_atom_begin(), r.related_atom_end() ); ! for( int i=0; i < ratoms.size(); ++i ) ! { ! remove_atom( ratoms[i] ); ! } ! remove(r); } *************** namespace mort *** 415,424 **** --- 458,478 ---- m_connecteds.clear(); } + void print_vec( std::ostream& os, const vector& v ) + { + for( int i=0; i < v.size(); ++i ) + { + os << v[i] << " "; + } + } + void molecule_t::cleanup( ) { map< hashid_t, map > idicts; + + // cleanup each component. meanwhile make a map between old ID and new ID, + // each component will have such a map, thus we use a map of map here. map< hashid_t, component_t >::iterator ci = m_components.begin(); for( ; ci != m_components.end(); ++ci ) { *************** namespace mort *** 437,452 **** --- 491,510 ---- assert( idicts.count(cid_1) > 0 && idicts.count(cid_2) > 0 ); + // idict_1 is the ID map of the component 1. idict_2 is the ID map of the component 2. map& idict_1 = idicts[cid_1]; map& idict_2 = idicts[cid_2]; adjacency_t anew; + // for each adjacency (ai->second), create a new adjacency, at the end use the + // new old replace the old one. int nobj = ai->second.size(); for( int i=0; i < nobj; ++i ) { vector* pvec = ai->second.getvec( i ); + // if idmap do not a ID, this mo must already be deleted, its relation should already be deleted. if( idict_1.count(i) == 0 ) { assert( pvec->size() == 0 ); Index: mortsrc/object/molecule.hpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/object/molecule.hpp,v retrieving revision 1.16 diff -p -r1.16 molecule.hpp *** src/gleap/mortsrc/object/molecule.hpp 30 Aug 2007 03:21:19 -0000 1.16 --- src/gleap/mortsrc/object/molecule.hpp 29 Apr 2008 19:20:55 -0000 *************** namespace mort *** 124,132 **** moref_t create_ptor(); moref_t create_oops(); void remove(const moref_t& mo); void remove_atom(const moref_t& a); ! void remove_resd(const moref_t& b); void clear(); void cleanup(); --- 124,136 ---- moref_t create_ptor(); moref_t create_oops(); + moref_t insert_atom( int relid ); + moref_t insert_resd( int relid ); + void remove(const moref_t& mo); void remove_atom(const moref_t& a); ! void remove_bond(const moref_t& b); ! void remove_resd(const moref_t& r); void clear(); void cleanup(); Index: mortsrc/object/molfun.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/object/molfun.cpp,v retrieving revision 1.14 diff -p -r1.14 molfun.cpp *** src/gleap/mortsrc/object/molfun.cpp 31 Aug 2007 14:41:27 -0000 1.14 --- src/gleap/mortsrc/object/molfun.cpp 29 Apr 2008 19:20:55 -0000 *************** namespace mort *** 37,43 **** return sum / m.natom(); } ! void minmaxpos( const molecule_t& m, numvec& pmin, numvec& pmax, bool vdwr ) { pmin.resize(3); pmax.resize(3); --- 37,43 ---- return sum / m.natom(); } ! void minmaxpos( const molecule_t& m, numvec& pmin, numvec& pmax, bool vdwr, int natom ) { pmin.resize(3); pmax.resize(3); *************** namespace mort *** 51,63 **** pmax[1] = pos[1]; pmax[2] = pos[2]; ! ! int natom = m.natom(); moiter_t ai = m.atom_begin(); for( int i=0; i < natom; ++i, ++ai ) { double add = vdwr ? ai->get_d(VDWR) : 0.0; - pmin[0] = std::min( pmin[0], pos[3*i ]-add ); pmin[1] = std::min( pmin[1], pos[3*i+1]-add ); pmin[2] = std::min( pmin[2], pos[3*i+2]-add ); --- 51,65 ---- pmax[1] = pos[1]; pmax[2] = pos[2]; ! if( natom==0 ) ! { ! natom = m.natom(); ! } ! moiter_t ai = m.atom_begin(); for( int i=0; i < natom; ++i, ++ai ) { double add = vdwr ? ai->get_d(VDWR) : 0.0; pmin[0] = std::min( pmin[0], pos[3*i ]-add ); pmin[1] = std::min( pmin[1], pos[3*i+1]-add ); pmin[2] = std::min( pmin[2], pos[3*i+2]-add ); *************** namespace mort *** 193,198 **** --- 195,206 ---- { const component_t* c = m.getcmp(cid); return c->get_vptr(parmid, 0); + } + + const string* get_sptr(const molecule_t& m, const hashid_t& cid, const hashid_t& parmid) + { + const component_t* c = m.getcmp(cid); + return c->get_sptr(parmid, 0); } matrix inter_moment( const molecule_t& m ) Index: mortsrc/object/molfun.hpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/object/molfun.hpp,v retrieving revision 1.10 diff -p -r1.10 molfun.hpp *** src/gleap/mortsrc/object/molfun.hpp 30 Aug 2007 03:20:46 -0000 1.10 --- src/gleap/mortsrc/object/molfun.hpp 29 Apr 2008 19:20:55 -0000 *************** namespace mort *** 17,23 **** numvec center( const molecule_t& m ); ! void minmaxpos( const molecule_t& m, numvec& pmin, numvec& pmax, bool vdwr=false ); void move( molecule_t& m, const numvec& v ); --- 17,23 ---- numvec center( const molecule_t& m ); ! void minmaxpos( const molecule_t& m, numvec& pmin, numvec& pmax, bool vdwr=false, int natom=0 ); void move( molecule_t& m, const numvec& v ); *************** namespace mort *** 30,35 **** --- 30,37 ---- const double* get_vptr(const molecule_t& m, const hashid_t& cid, const hashid_t& parmid); const double* get_dptr(const molecule_t& m, const hashid_t& cid, const hashid_t& parmid); + + const string* get_sptr(const molecule_t& m, const hashid_t& cid, const hashid_t& parmid); const int* get_iptr(const molecule_t& m, const hashid_t& cid, const hashid_t& parmid); Index: mortsrc/object/resdfun.cpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/object/resdfun.cpp,v retrieving revision 1.14 diff -p -r1.14 resdfun.cpp *** src/gleap/mortsrc/object/resdfun.cpp 1 Jul 2007 16:27:54 -0000 1.14 --- src/gleap/mortsrc/object/resdfun.cpp 29 Apr 2008 19:20:55 -0000 *************** namespace mort *** 45,67 **** return sum; } ! moref_t merge( molecule_t& mol, const molecule_t& sub) { ! moref_t resd = mol.create_resd(); ! ! moiter_t a = sub.atom_begin(); ! for( ; a != sub.atom_end(); ++a ) ! { ! build_atom(resd, *a); ! } ! ! ! moiter_t b = sub.bond_begin(); ! for( ; b != sub.bond_end(); ++b ) ! { ! build_bond( resd, *b); ! } ! string name; if( sub.get_s(NAME, name) ) --- 45,53 ---- return sum; } ! moref_t merge( molecule_t& mol, const molecule_t& sub, int relid) { ! moref_t resd = (relid<0) ? mol.create_resd() : mol.insert_resd(relid); string name; if( sub.get_s(NAME, name) ) *************** namespace mort *** 75,80 **** --- 61,78 ---- resd.set_s(TYPE, type); } + moiter_t a = sub.atom_begin(); + for( ; a != sub.atom_end(); ++a ) + { + build_atom(resd, *a); + } + + moiter_t b = sub.bond_begin(); + for( ; b != sub.bond_end(); ++b ) + { + build_bond( resd, *b); + } + int head; if( sub.get_i(HEAD,head) ) { *************** namespace mort *** 137,143 **** return resd; } ! moref_t build_atom(moref_t& target, const moref_t& atomsrc, const namemap_t* nmap) { vector excepts; numvec oldpos; --- 135,141 ---- return resd; } ! moref_t build_atom(moref_t& resddst, const moref_t& atomsrc, const namemap_t* nmap) { vector excepts; numvec oldpos; *************** namespace mort *** 148,194 **** atom_name = nmap->get_atom_name(atom_name); } ! moref_t dst( target.getmol(), ATOM, -1); ! if( get_atom(target, atom_name, dst) ) { ! dst.set_i(TYPEID, KNOWN); ! dst.get_v(POSITION, oldpos); } else { ! dst = create_atom(target, atom_name); ! dst.set_i(TYPEID, UNKNOWN); } ! copy_allparm_except(atomsrc, dst, excepts); ! if( dst.get_i(TYPEID)==(int)KNOWN ) { ! dst.set_v(POSITION, oldpos); } ! return dst; } ! moref_t build_bond(moref_t& target, const moref_t& bondsrc) { string name_0 = atom_1st( bondsrc ).get_s(NAME); string name_1 = atom_2nd( bondsrc ).get_s(NAME); ! moref_t a_0 = get_atom(target, name_0 ); ! moref_t a_1 = get_atom(target, name_1 ); ! moref_t out = frcget_bond( a_0, a_1 ); ! copy_allparm(bondsrc, out); ! return out; } bool get_atom(const moref_t& resd, const string& name, moref_t& atom) { return resd.get_related_atom( sparm_comper1(NAME, name), atom); } moref_t get_atom(const moref_t& resd, const string& name) { return resd.get_related_atom( sparm_comper1(NAME, name) ); } --- 146,194 ---- atom_name = nmap->get_atom_name(atom_name); } ! moref_t atomdst( resddst.getmol(), ATOM, -1); ! if( get_atom(resddst, atom_name, atomdst) ) { ! atomdst.set_i(TYPEID, KNOWN); ! atomdst.get_v(POSITION, oldpos); } else { ! atomdst = create_atom(resddst, atom_name); ! atomdst.set_i(TYPEID, UNKNOWN); } ! copy_allparm_except(atomsrc, atomdst, excepts); ! if( atomdst.get_i(TYPEID)==(int)KNOWN ) { ! atomdst.set_v(POSITION, oldpos); } ! return atomdst; } ! moref_t build_bond(moref_t& resddst, const moref_t& bondsrc) { string name_0 = atom_1st( bondsrc ).get_s(NAME); string name_1 = atom_2nd( bondsrc ).get_s(NAME); ! moref_t a_0 = get_atom(resddst, name_0 ); ! moref_t a_1 = get_atom(resddst, name_1 ); ! moref_t bonddst = frcget_bond( a_0, a_1 ); ! copy_allparm(bondsrc, bonddst); ! return bonddst; } bool get_atom(const moref_t& resd, const string& name, moref_t& atom) { + assert( resd.cmpid()==RESD ); return resd.get_related_atom( sparm_comper1(NAME, name), atom); } moref_t get_atom(const moref_t& resd, const string& name) { + assert( resd.cmpid()==RESD ); return resd.get_related_atom( sparm_comper1(NAME, name) ); } *************** namespace mort *** 198,209 **** return get_atom(resd, name, atom) ? atom : create_atom( resd, name ); } moref_t create_atom(moref_t& r, const string& name) { ! moref_t a = r.getmol().create_atom(); a.connect( r ); r.connect( a ); - a.set_s(NAME, name); return a; } --- 198,234 ---- return get_atom(resd, name, atom) ? atom : create_atom( resd, name ); } + bool is_last_resd( moref_t& r ) + { + assert( r.cmpid()==RESD ); + return r.relid()==r.getmol().nresd()-1; + } + + // for last residue: insert point is natom, + // otherwise it is the relative ID of the first atom of the next residue. + int insert_point( moref_t& r ) + { + if( is_last_resd(r) ) + { + return r.getmol().natom(); + } + + molecule_t* pmol = &r.getmol(); + int resdpos = r.relid(); + moiter_t rnext = pmol->resd_begin() + resdpos + 1; + assert( rnext != pmol->resd_end() ); + + return rnext->related_atom_begin()->relid(); + } + moref_t create_atom(moref_t& r, const string& name) { ! molecule_t* pmol = &r.getmol(); ! moref_t a = is_last_resd(r) ? pmol->create_atom() : pmol->insert_atom( insert_point(r) ); ! a.set_s(NAME, name); ! a.connect( r ); r.connect( a ); return a; } *************** namespace mort *** 265,270 **** --- 290,312 ---- } + bool has_interbond( const moref_t& ri, const moref_t& rj ) + { + assert( ri.cmpid()==RESD && rj.cmpid()==RESD ); + + moiter_t ai = ri.related_atom_begin(); + for( ; ai != ri.related_atom_end(); ++ai ) + { + moiter_t aj = ai->related_atom_begin(); + for( ; aj != ai->related_atom_end(); ++aj ) + { + if( aj->get_i(RESID)==rj.get_i(ID) ) + return true; + } + } + + return false; + } } Index: mortsrc/object/resdfun.hpp =================================================================== RCS file: /raid5/case/cvsroot/amber10/src/gleap/mortsrc/object/resdfun.hpp,v retrieving revision 1.6 diff -p -r1.6 resdfun.hpp *** src/gleap/mortsrc/object/resdfun.hpp 5 May 2007 20:27:54 -0000 1.6 --- src/gleap/mortsrc/object/resdfun.hpp 29 Apr 2008 19:20:55 -0000 *************** namespace mort *** 21,27 **** moref_t merge( molecule_t& mol, const moref_t& subresd ); ! moref_t merge( molecule_t& mol, const molecule_t& submol ); moref_t build_atom(moref_t& resd, const moref_t& atomsrc, const namemap_t* nmap=NULL); --- 21,27 ---- moref_t merge( molecule_t& mol, const moref_t& subresd ); ! moref_t merge( molecule_t& mol, const molecule_t& submol, int relid=-1 ); moref_t build_atom(moref_t& resd, const moref_t& atomsrc, const namemap_t* nmap=NULL); *************** namespace mort *** 39,44 **** --- 39,45 ---- moref_t copy_resd( const moref_t& r, molecule_t& m ); + bool has_interbond( const moref_t& ri, const moref_t& rj ); }; #endif --------------------------------------------------------------------------- Workarounds: none ********>Bugfix 7 Author: Jianyin Shao Date: 06/05/2008 Programs: ptraj Description: The bug in PtrajClusteringGetAttributeArrayCount() of cluster.c may point to an empty pointer, therefore could cause a runtime error. The other bug in mask.c will falsely issue a warning for "*" for empty mask, this bug will not cause any error. Fix: This patch acctually affects multiple files. So, run it from the $AMBERHOME directory, with: $ patch -p0 -N -r patch_rejects < bugfix.6 ------------------------------------------------------------------------------ diff -Naur oldamber/src/ptraj/cluster.c newamber/src/ptraj/cluster.c --- src/ptraj/cluster.c 2008-02-22 10:52:34.000000000 -0700 +++ src/ptraj/cluster.c 2008-06-05 15:13:21.000000000 -0600 @@ -6513,7 +6513,7 @@ Cluster* NewCluster; int PointInCluster; int ClusterCount = This->ClusterCount; - int AttributeCount = PtrajClusteringGetAttributeArrayCount(This); + int AttributeCount = PtrajClusteringGetAttributeCount(This); int i, j, k, diff; int* AttributeDifference; @@ -6806,8 +6806,8 @@ ClusterFindBestP2C( (PtrajCluster *) Node->Cluster); Node = Node->pNext; } - PtrajClusteringSOMExtractDifference( (PtrajClustering *) This, SOMNodes); /* print out significant differences in attributes between clusters. */ + if (prnlev >6) PtrajClusteringSOMExtractDifference( (PtrajClustering *) This, SOMNodes); /* print out significant differences in attributes between clusters. */ } /* diff -Naur oldamber/src/ptraj/cluster.c newamber/src/ptraj/mask.c --- src/ptraj/mask.c 2008-02-22 10:52:34.000000000 -0700 +++ src/ptraj/mask.c 2008-06-05 15:13:21.000000000 -0600 @@ -397,8 +397,8 @@ flag = 1; n++; } else if ( *p == ']' ) { - if ( n < 3 ) - printf("Warning: empty token?\n"); + if ( n < 3 && *(p-1) != '*' ) + printf("Warning: \'%c\' empty token?\n", *(p-1)); flag = 0; n = 0; } else { ------------------------------------------------------------------------------ Temporary workarounds: none