From 3dd30109ec4c63e67949ba74e7e057264df65529 Mon Sep 17 00:00:00 2001 From: Gabriel Studer <gabriel.studer@unibas.ch> Date: Mon, 10 Aug 2015 18:06:41 +0200 Subject: [PATCH] add hydrogen bond term as described in FG-MD paper from Zhang --- modules/mol/mm/pymod/export_topology.cc | 46 +++++++++ modules/mol/mm/src/system_creator.cc | 60 +++++++++++ modules/mol/mm/src/topology.cc | 132 ++++++++++++++++++++++++ modules/mol/mm/src/topology.hh | 69 ++++++++++++- 4 files changed, 305 insertions(+), 2 deletions(-) diff --git a/modules/mol/mm/pymod/export_topology.cc b/modules/mol/mm/pymod/export_topology.cc index 9b26eed6a..f72bd7502 100644 --- a/modules/mol/mm/pymod/export_topology.cc +++ b/modules/mol/mm/pymod/export_topology.cc @@ -218,6 +218,19 @@ namespace{ return boost::python::make_tuple(i1,i2,l,k); } + boost::python::tuple WrapGetFGMDHBondDonorParam(ost::mol::mm::TopologyPtr p, uint index){ + uint i1, i2; + Real d,k_d,a,k_a,b,k_b; + p->GetFGMDHBondDonorParameters(index,i1,i2,d,k_d,a,k_a,b,k_b); + return boost::python::make_tuple(i1,i2,d,k_d,a,k_a,b,k_b); + } + + boost::python::tuple WrapGetFGMDHBondAcceptorParam(ost::mol::mm::TopologyPtr p, uint index){ + uint i1, i2; + p->GetFGMDHBondAcceptorParameters(index,i1,i2); + return boost::python::make_tuple(i1,i2); + } + void WrapSetCMapParameters(ost::mol::mm::TopologyPtr p, uint index, int dimension, boost::python::list l){ std::vector<Real> v = ListToVec<Real>(l); p->SetCMapParameters(index,dimension,v); @@ -313,11 +326,32 @@ namespace{ return VecToList<uint>(return_vec); } + boost::python::list GetFGMDHBondDonorIndices(ost::mol::mm::TopologyPtr p, uint i1, uint i2){ + std::vector<uint> return_vec = p->GetFGMDHBondDonorIndices(i1,i2); + return VecToList<uint>(return_vec); + } + + boost::python::list GetFGMDHBondDonorIndicesSingleIndex(ost::mol::mm::TopologyPtr p, uint i){ + std::vector<uint> return_vec = p->GetFGMDHBondDonorIndices(i); + return VecToList<uint>(return_vec); + } + boost::python::list GetHarmonicPositionRestraintIndicesSingleIndex(ost::mol::mm::TopologyPtr p, uint i){ std::vector<uint> return_vec = p->GetHarmonicPositionRestraintIndices(i); return VecToList<uint>(return_vec); } + boost::python::list GetFGMDHBondAcceptorIndices(ost::mol::mm::TopologyPtr p, uint i1, uint i2){ + std::vector<uint> return_vec = p->GetFGMDHBondAcceptorIndices(i1,i2); + return VecToList<uint>(return_vec); + } + + boost::python::list GetFGMDHBondAcceptorIndicesSingleIndex(ost::mol::mm::TopologyPtr p, uint i){ + std::vector<uint> return_vec = p->GetFGMDHBondAcceptorIndices(i); + return VecToList<uint>(return_vec); + } + + void MergeTop(ost::mol::mm::TopologyPtr top, ost::mol::mm::TopologyPtr other){ top->Merge(other); } @@ -357,6 +391,8 @@ void export_Topology() .def("ResetExclusions",&ost::mol::mm::Topology::ResetExclusions) .def("AddHarmonicPositionRestraint",&ost::mol::mm::Topology::AddHarmonicPositionRestraint,(arg("index"),arg("ref_position"),arg("k"),arg("x_scale")=1.0,arg("y_scale")=1.0,arg("z_scale")=1.0)) .def("AddHarmonicDistanceRestraint",&ost::mol::mm::Topology::AddHarmonicDistanceRestraint) + .def("AddFGMDHBondDonor",&ost::mol::mm::Topology::AddFGMDHBondDonor) + .def("AddFGMDHBondAcceptor",&ost::mol::mm::Topology::AddFGMDHBondAcceptor) //single atom parameter getter and setter functions .def("SetSigmas",&WrapSetSigmas) @@ -400,6 +436,8 @@ void export_Topology() .def("GetDistanceConstraintParameters",&WrapGetDistanceConstraintParam) .def("GetHarmonicPositionRestraintParameters",&WrapGetHarmonicPositionRestraintParam) .def("GetHarmonicDistanceRestraintParameters",&WrapGetHarmonicDistanceRestraintParam) + .def("GetFGMDHBondDonorParameters",&WrapGetFGMDHBondDonorParam) + .def("GetFGMDHBondAcceptorParameters",&WrapGetFGMDHBondAcceptorParam) //setter functions for interaction parameters .def("SetHarmonicBondParameters",&ost::mol::mm::Topology::SetHarmonicBondParameters) @@ -413,6 +451,7 @@ void export_Topology() .def("SetDistanceConstraintParameters",&ost::mol::mm::Topology::SetDistanceConstraintParameters) .def("SetHarmonicPositionRestraintParameters",&ost::mol::mm::Topology::SetHarmonicPositionRestraintParameters) .def("SetHarmonicDistanceRestraintParameters",&ost::mol::mm::Topology::SetHarmonicDistanceRestraintParameters) + .def("SetFGMDHBondDonorParameters",&ost::mol::mm::Topology::SetFGMDHBondDonorParameters) //functions to find interactions certain atoms are involved in .def("GetHarmonicBondIndices",&GetHarmonicBondIndices) @@ -436,6 +475,11 @@ void export_Topology() .def("GetHarmonicDistanceRestraintIndices",&GetHarmonicDistanceRestraintIndices) .def("GetHarmonicDistanceRestraintIndices",&GetHarmonicDistanceRestraintIndicesSingleIndex) .def("GetHarmonicPositionRestraintIndices",&GetHarmonicPositionRestraintIndicesSingleIndex) + .def("GetFGMDHBondDonorIndices",&GetFGMDHBondDonorIndices) + .def("GetFGMDHBondDonorIndices",&GetFGMDHBondDonorIndicesSingleIndex) + .def("GetFGMDHBondAcceptorIndices",&GetFGMDHBondAcceptorIndices) + .def("GetFGMDHBondAcceptorIndices",&GetFGMDHBondAcceptorIndicesSingleIndex) + //functions to get amount of data in topology .def("GetNumParticles",&ost::mol::mm::Topology::GetNumParticles) @@ -450,6 +494,8 @@ void export_Topology() .def("GetNumExclusions",&ost::mol::mm::Topology::GetNumExclusions) .def("GetNumHarmonicPositionRestraints",&ost::mol::mm::Topology::GetNumHarmonicPositionRestraints) .def("GetNumHarmonicDistanceRestraints",&ost::mol::mm::Topology::GetNumHarmonicDistanceRestraints) + .def("GetNumFGMDHBondDonors",&ost::mol::mm::Topology::GetNumFGMDHBondDonors) + .def("GetNumFGMDHBondAcceptors",&ost::mol::mm::Topology::GetNumFGMDHBondAcceptors) .def("Merge",&MergeTop) .def("Merge",&MergeTopEnt) diff --git a/modules/mol/mm/src/system_creator.cc b/modules/mol/mm/src/system_creator.cc index 5394ddf86..2912a5dca 100644 --- a/modules/mol/mm/src/system_creator.cc +++ b/modules/mol/mm/src/system_creator.cc @@ -236,6 +236,66 @@ SystemPtr SystemCreator::Create(const TopologyPtr top, } } + const std::vector<std::pair<Index<2>,std::vector<Real> > > fgmd_hbond_donors = top->GetFGMDHBondDonors(); + const std::vector<Index<2> > fgmd_hbond_acceptors = top->GetFGMDHBondAcceptors(); + if(!fgmd_hbond_donors.empty() && !fgmd_hbond_acceptors.empty()){ + String distance_term = "k_length*(distance(a1,d1)-k_length)^2"; + String alpha_term = "k_alpha*(angle(d2,d1,a1)-alpha)^2"; + String beta_term = "k_beta*(angle(d1,a1,a2)-beta)^2"; + String f = distance_term + "+" + alpha_term + "+" + beta_term; + OpenMM::CustomHbondForce& fgmd_hbond_force = *new OpenMM::CustomHbondForce(f); + + fgmd_hbond_force.addPerDonorParameter("length"); + fgmd_hbond_force.addPerDonorParameter("k_length"); + fgmd_hbond_force.addPerDonorParameter("alpha"); + fgmd_hbond_force.addPerDonorParameter("k_alpha"); + fgmd_hbond_force.addPerDonorParameter("beta"); + fgmd_hbond_force.addPerDonorParameter("k_beta"); + + fgmd_hbond_force.setCutoffDistance(3.0); + fgmd_hbond_force.setNonbondedMethod(OpenMM::CustomHbondForce::CutoffPeriodic); + + switch(settings->nonbonded_method){ + case NoCutoff:{ + fgmd_hbond_force.setNonbondedMethod(OpenMM::CustomHbondForce::NoCutoff); + break; + } + case CutoffNonPeriodic:{ + fgmd_hbond_force.setNonbondedMethod(OpenMM::CustomHbondForce::CutoffNonPeriodic); + break; + } + case CutoffPeriodic:{ + fgmd_hbond_force.setNonbondedMethod(OpenMM::CustomHbondForce::CutoffPeriodic); + break; + } + case Ewald:{ + fgmd_hbond_force.setNonbondedMethod(OpenMM::CustomHbondForce::CutoffPeriodic); + break; + } + case PME:{ + fgmd_hbond_force.setNonbondedMethod(OpenMM::CustomHbondForce::CutoffPeriodic); + break; + } + } + + for(std::vector<std::pair<Index<2>,std::vector<Real> > >::const_iterator i = fgmd_hbond_donors.begin(); + i != fgmd_hbond_donors.end(); ++i){ + //stupid cast + std::vector<double> parameters; + for(std::vector<Real>::const_iterator j = i->second.begin(); + j != i->second.end(); ++j){ + parameters.push_back(*j); + } + fgmd_hbond_force.addDonor(i->first[0],i->first[1],-1, parameters); + } + std::vector<double> dummy_vec; + for(std::vector<Index<2> >::const_iterator i = fgmd_hbond_acceptors.begin(); + i != fgmd_hbond_acceptors.end(); ++i){ + fgmd_hbond_force.addAcceptor((*i)[0],(*i)[1],-1,dummy_vec); + } + sys->addForce(&fgmd_hbond_force); + } + std::vector<Real> sigmas = top->GetSigmas(); std::vector<Real> epsilons = top->GetEpsilons(); std::vector<Real> charges = top->GetCharges(); diff --git a/modules/mol/mm/src/topology.cc b/modules/mol/mm/src/topology.cc index 3a3676449..448e25490 100644 --- a/modules/mol/mm/src/topology.cc +++ b/modules/mol/mm/src/topology.cc @@ -287,6 +287,33 @@ uint Topology::AddHarmonicDistanceRestraint(uint index_one, uint index_two, return harmonic_distance_restraints_.size()-1; } +uint Topology::AddFGMDHBondDonor(uint index_one, uint index_two, + Real length, Real k_length, Real alpha, + Real k_alpha, Real beta, Real k_beta){ + if(index_one >= num_particles_ || index_two >= num_particles_){ + throw ost::Error("Index of fgmd hbond donor atom exceeds number of particles present in topology!"); + } + Index<2> index(index_one,index_two); + std::vector<Real> parameters; + parameters.push_back(length); + parameters.push_back(k_length); + parameters.push_back(alpha); + parameters.push_back(k_alpha); + parameters.push_back(beta); + parameters.push_back(k_beta); + fgmd_hbond_donors_.push_back(std::make_pair(index,parameters)); + return fgmd_hbond_donors_.size()-1; +} + +uint Topology::AddFGMDHBondAcceptor(uint index_one, uint index_two){ + if(index_one >= num_particles_ || index_two >= num_particles_){ + throw ost::Error("Index of fgmd hbond acceptor atom exceeds number of particles present in topology!"); + } + Index<2> index(index_one,index_two); + fgmd_hbond_acceptors_.push_back(index); + return fgmd_hbond_acceptors_.size()-1; +} + void Topology::SetSigmas(const std::vector<Real>& sigmas){ if(sigmas.size() != num_particles_){ throw ost::Error("Expect the same number of sigma parameters than particles!"); @@ -531,6 +558,29 @@ void Topology::GetHarmonicDistanceRestraintParameters(uint index, uint& atom_one force_constant = harmonic_distance_restraints_[index].second[1]; } +void Topology::GetFGMDHBondDonorParameters(uint index, uint& atom_one, uint& atom_two, Real& length, + Real& k_length, Real& alpha, Real& k_alpha, + Real& beta, Real& k_beta) const{ + if(index >= fgmd_hbond_donors_.size()){ + throw ost::Error("Provided index exceeds number of fgmd hbond donors present in topology!"); + } + atom_one = harmonic_distance_restraints_[index].first[0]; + atom_two = harmonic_distance_restraints_[index].first[1]; + length = fgmd_hbond_donors_[index].second[0]; + k_length = fgmd_hbond_donors_[index].second[1]; + alpha = fgmd_hbond_donors_[index].second[2]; + k_alpha = fgmd_hbond_donors_[index].second[3]; + beta = fgmd_hbond_donors_[index].second[4]; + k_beta = fgmd_hbond_donors_[index].second[5]; +} + +void Topology::GetFGMDHBondAcceptorParameters(uint index, uint& atom_one, uint& atom_two) const{ + if(index >= fgmd_hbond_acceptors_.size()){ + throw ost::Error("Provided index exceeds number of fgmd hbond acceptors present in topology!"); + } + atom_one = fgmd_hbond_acceptors_[index][0]; + atom_two = fgmd_hbond_acceptors_[index][1]; +} void Topology::SetHarmonicBondParameters(uint index, const Real bond_length, const Real force_constant){ if(index >= harmonic_bonds_.size()){ @@ -630,7 +680,19 @@ void Topology::SetHarmonicDistanceRestraintParameters(uint index, Real length, } harmonic_distance_restraints_[index].second[0] = length; harmonic_distance_restraints_[index].second[1] = force_constant; +} +void Topology::SetFGMDHBondDonorParameters(uint index, Real length, Real k_length, + Real alpha, Real k_alpha, Real beta, Real k_beta){ + if(index >= fgmd_hbond_donors_.size()){ + throw ost::Error("Provided index exceeds number of fgmd donors present in topology!"); + } + fgmd_hbond_donors_[index].second[0] = length; + fgmd_hbond_donors_[index].second[1] = k_length; + fgmd_hbond_donors_[index].second[2] = alpha; + fgmd_hbond_donors_[index].second[3] = k_alpha; + fgmd_hbond_donors_[index].second[4] = beta; + fgmd_hbond_donors_[index].second[5] = k_beta; } Real Topology::GetMass(uint index) const{ @@ -819,6 +881,32 @@ std::vector<uint> Topology::GetHarmonicDistanceRestraintIndices(uint index_one, return return_indices; } +std::vector<uint> Topology::GetFGMDHBondDonorIndices(uint index_one, uint index_two) const{ + + Index<2> index(index_one,index_two); + Index<2> reverse_index(index_two,index_one); + std::vector<uint> return_indices; + for(uint i = 0; i < fgmd_hbond_donors_.size(); ++i){ + if(index == fgmd_hbond_donors_[i].first || reverse_index == fgmd_hbond_donors_[i].first){ + return_indices.push_back(i); + } + } + return return_indices; +} + +std::vector<uint> Topology::GetFGMDHBondAcceptorIndices(uint index_one, uint index_two) const{ + + Index<2> index(index_one,index_two); + Index<2> reverse_index(index_two,index_one); + std::vector<uint> return_indices; + for(uint i = 0; i < fgmd_hbond_acceptors_.size(); ++i){ + if(index == fgmd_hbond_acceptors_[i] || reverse_index == fgmd_hbond_acceptors_[i]){ + return_indices.push_back(i); + } + } + return return_indices; +} + std::vector<uint> Topology::GetHarmonicBondIndices(uint atom_index) const{ std::vector<uint> return_indices; for(uint i = 0; i < harmonic_bonds_.size(); ++i){ @@ -926,6 +1014,26 @@ std::vector<uint> Topology::GetHarmonicDistanceRestraintIndices(uint atom_index) return return_indices; } +std::vector<uint> Topology::GetFGMDHBondDonorIndices(uint atom_index) const{ + std::vector<uint> return_indices; + for(uint i = 0; i < fgmd_hbond_donors_.size(); ++i){ + if(fgmd_hbond_donors_[i].first[0] == atom_index || fgmd_hbond_donors_[i].first[1] == atom_index){ + return_indices.push_back(i); + } + } + return return_indices; +} + +std::vector<uint> Topology::GetFGMDHBondAcceptorIndices(uint atom_index) const{ + std::vector<uint> return_indices; + for(uint i = 0; i < fgmd_hbond_acceptors_.size(); ++i){ + if(fgmd_hbond_acceptors_[i][0] == atom_index || fgmd_hbond_acceptors_[i][1] == atom_index){ + return_indices.push_back(i); + } + } + return return_indices; +} + void Topology::Merge(ost::mol::EntityHandle& ent, TopologyPtr other, const ost::mol::EntityHandle& other_ent){ @@ -1116,6 +1224,30 @@ void Topology::MergeTop(TopologyPtr other){ } } + const std::vector<std::pair<Index<2>,std::vector<Real> > >& donors = other->GetFGMDHBondDonors(); + if(!donors.empty()){ + for(std::vector<std::pair<Index<2>,std::vector<Real> > >::const_iterator i = donors.begin(); + i != donors.end(); ++i){ + this->AddFGMDHBondDonor(old_num_particles + i->first[0], + old_num_particles + i->first[1], + i->second[0], + i->second[1], + i->second[2], + i->second[3], + i->second[4], + i->second[5]); + } + } + + const std::vector<Index<2> >& acceptors = other->GetFGMDHBondAcceptors(); + if(!acceptors.empty()){ + for(std::vector<Index<2> >::const_iterator i = acceptors.begin(); + i != acceptors.end(); ++i){ + this->AddFGMDHBondAcceptor(old_num_particles + (*i)[0], + old_num_particles + (*i)[1]); + } + } + const std::vector<uint>& position_constraints = other->GetPositionConstraints(); if(!position_constraints.empty()){ for(std::vector<uint>::const_iterator i = position_constraints.begin(); diff --git a/modules/mol/mm/src/topology.hh b/modules/mol/mm/src/topology.hh index b626169fa..e578e4200 100644 --- a/modules/mol/mm/src/topology.hh +++ b/modules/mol/mm/src/topology.hh @@ -137,6 +137,12 @@ public: uint AddHarmonicDistanceRestraint(uint index_one, uint index_two, Real length, Real force_constant); + uint AddFGMDHBondDonor(uint index_one, uint index_two, + Real length, Real k_length, Real alpha, + Real k_alpha, Real beta, Real k_beta); + + uint AddFGMDHBondAcceptor(uint index_one, uint index_two); + //Single atom parameters are expected to be set at once... void SetSigmas(const std::vector<Real>& sigmas); @@ -195,10 +201,16 @@ public: Real& distance) const; void GetHarmonicPositionRestraintParameters(uint index, uint& atom_index, geom::Vec3& ref_position, - Real& k, Real& x_scale, Real& y_scale, Real& z_scale) const; + Real& k, Real& x_scale, Real& y_scale, Real& z_scale) const; void GetHarmonicDistanceRestraintParameters(uint index, uint& atom_one, uint& atom_two, Real& length, - Real& force_constant) const; + Real& force_constant) const; + + void GetFGMDHBondDonorParameters(uint index, uint& index_one, uint& index_two, + Real& length, Real& k_length, Real& alpha, + Real& k_alpha, Real& beta, Real& k_beta) const; + + void GetFGMDHBondAcceptorParameters(uint index, uint& index_one, uint& index_two) const; void SetHarmonicBondParameters(uint index, const Real bond_length, const Real force_constant); @@ -226,6 +238,10 @@ public: void SetHarmonicDistanceRestraintParameters(uint index, Real length, Real force_constant); + void SetFGMDHBondDonorParameters(uint index, Real length, Real k_length, + Real alpha, Real k_alpha, Real beta, + Real k_beta); + const std::vector<std::pair<Index<2>, std::vector<Real> > >& GetHarmonicBonds() const { return harmonic_bonds_; } const std::vector<std::pair<Index<3>, std::vector<Real> > >& GetHarmonicAngles() const { return harmonic_angles_; } @@ -252,6 +268,10 @@ public: const std::vector<uint>& GetPositionConstraints() const { return position_constraints_; } + const std::vector<std::pair<Index<2>, std::vector<Real> > >& GetFGMDHBondDonors() const { return fgmd_hbond_donors_; } + + const std::vector<Index<2> >& GetFGMDHBondAcceptors() const { return fgmd_hbond_acceptors_; } + std::vector<Real> GetSigmas() const { return sigmas_; } std::vector<Real> GetEpsilons() const { return epsilons_; } @@ -323,6 +343,11 @@ public: std::vector<uint> GetHarmonicDistanceRestraintIndices(uint index_one, uint index_two) const; + std::vector<uint> GetFGMDHBondDonorIndices(uint index_one, + uint index_two) const; + + std::vector<uint> GetFGMDHBondAcceptorIndices(uint index_one, + uint index_two) const; std::vector<uint> GetHarmonicBondIndices(uint atom_index) const; @@ -346,6 +371,10 @@ public: std::vector<uint> GetHarmonicDistanceRestraintIndices(uint atom_index) const; + std::vector<uint> GetFGMDHBondDonorIndices(uint atom_index) const; + + std::vector<uint> GetFGMDHBondAcceptorIndices(uint atom_index) const; + uint GetNumParticles() { return num_particles_; } uint GetNumHarmonicBonds() { return harmonic_bonds_.size(); } @@ -374,6 +403,10 @@ public: uint GetNumExclusions() { return exclusions_.size(); } + uint GetNumFGMDHBondDonors() { return fgmd_hbond_donors_.size(); } + + uint GetNumFGMDHBondAcceptors() { return fgmd_hbond_acceptors_.size(); } + void Merge(ost::mol::EntityHandle& ent, TopologyPtr other, const ost::mol::EntityHandle& other_ent); void Merge(TopologyPtr other); @@ -517,6 +550,16 @@ public: for(uint i = 0; i < num_items; ++i){ harmonic_distance_restraints_.push_back(std::make_pair(Index<2>(),std::vector<Real>(2))); } + + ds & num_items; + for(uint i = 0; i < num_items; ++i){ + fgmd_hbond_donors_.push_back(std::make_pair(Index<2>(),std::vector<Real>(6))); + } + + ds & num_items; + for(uint i = 0; i < num_items; ++i){ + fgmd_hbond_acceptors_.push_back(Index<2>()); + } } else{ num_items = harmonic_bonds_.size(); @@ -547,6 +590,10 @@ public: ds & num_items; num_items = harmonic_distance_restraints_.size(); ds & num_items; + num_items = fgmd_hbond_donors_.size(); + ds & num_items; + num_items = fgmd_hbond_acceptors_.size(); + ds & num_items; } for(std::vector<std::pair<Index<2>,std::vector<Real> > >::iterator i = harmonic_bonds_.begin(); @@ -641,6 +688,22 @@ public: ds & i->second[1]; } + for(std::vector<std::pair<Index<2>,std::vector<Real> > >::iterator i = fgmd_hbond_donors_.begin(); + i != fgmd_hbond_donors_.end(); ++i){ + ds & i->first; + ds & i->second[0]; + ds & i->second[1]; + ds & i->second[2]; + ds & i->second[3]; + ds & i->second[4]; + ds & i->second[5]; + } + + for(std::vector<Index<2> >::iterator i = fgmd_hbond_acceptors_.begin(); + i != fgmd_hbond_acceptors_.end(); ++i){ + ds & (*i); + } + if(ds.IsSource()){ ds & num_items; for(uint i = 0; i < num_items; ++i){ @@ -724,6 +787,8 @@ private: std::vector<Index<2> > exclusions_; std::vector<std::pair<Index<1>,std::vector<Real> > > harmonic_position_restraints_; std::vector<std::pair<Index<2>,std::vector<Real> > > harmonic_distance_restraints_; + std::vector<std::pair<Index<2>,std::vector<Real> > > fgmd_hbond_donors_; + std::vector<Index<2> > fgmd_hbond_acceptors_; //the atoms of the interactions, that should be unique get tracked in here //note, that this is waste of memory, needs better implementation -- GitLab