********> 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
*** mortsrc/capbox/addions.cpp	26 Jun 2007 13:44:00 -0000	1.10
--- 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
*** mortsrc/capbox/octree.cpp	21 Dec 2006 17:27:01 -0000	1.5
--- 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
*** mortsrc/capbox/solvate.cpp	9 Apr 2008 23:34:00 -0000	1.17
--- 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
*** mortsrc/capbox/solvate.hpp	9 Apr 2008 23:34:00 -0000	1.5
--- 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
*** mortsrc/common/hash.txt	21 Jan 2008 22:20:10 -0000	1.6
--- 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
*** mortsrc/common/hashcode.hpp	21 Jan 2008 22:20:10 -0000	1.10
--- 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
*** mortsrc/guilib/mainwin.cpp	2 Apr 2008 18:25:11 -0000	1.11
--- 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
*** mortsrc/object/component.cpp	31 Aug 2007 21:44:50 -0000	1.15
--- 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
*** mortsrc/object/component.hpp	4 May 2007 12:18:04 -0000	1.8
--- 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
*** mortsrc/object/database.cpp	10 Feb 2008 15:04:43 -0000	1.12
--- 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
*** mortsrc/object/database.hpp	17 Aug 2007 03:04:43 -0000	1.14
--- 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
*** mortsrc/object/entity.cpp	30 Aug 2007 03:23:02 -0000	1.9
--- 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
*** mortsrc/object/entity.hpp	30 Aug 2007 03:22:49 -0000	1.6
--- 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
*** mortsrc/object/molecule.cpp	20 Aug 2007 18:22:40 -0000	1.15
--- 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
*** mortsrc/object/molecule.hpp	30 Aug 2007 03:21:19 -0000	1.16
--- 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
*** mortsrc/object/molfun.cpp	31 Aug 2007 14:41:27 -0000	1.14
--- 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
*** mortsrc/object/molfun.hpp	30 Aug 2007 03:20:46 -0000	1.10
--- 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
*** mortsrc/object/resdfun.cpp	1 Jul 2007 16:27:54 -0000	1.14
--- 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
*** mortsrc/object/resdfun.hpp	5 May 2007 20:27:54 -0000	1.6
--- 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