********> 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.
*** mortsrc/capbox.hpp 18 Dec 2006 14:54:52 -0000 1.4
--- 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
*** mortsrc/capbox/solvate.cpp 21 Jan 2008 22:21:22 -0000 1.16
--- 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
*** mortsrc/capbox/solvate.hpp 21 Dec 2006 17:27:15 -0000 1.4
--- 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
*** plugins/solvate.cpp 9 May 2007 19:46:06 -0000 1.7
--- 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
*** plugins/solvate.hpp 9 May 2007 19:46:06 -0000 1.2
--- 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