Something went wrong on our end
-
Gerardo Tauriello authoredGerardo Tauriello authored
gromacs_block_modifiers.cc 14.95 KiB
//------------------------------------------------------------------------------
// This file is part of the OpenStructure project <www.openstructure.org>
//
// Copyright (C) 2008-2015 by the OpenStructure authors
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License as published by the Free
// Software Foundation; either version 3.0 of the License, or (at your option)
// any later version.
// This library is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
// details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this library; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
//------------------------------------------------------------------------------
#include <ost/mol/mm/gromacs_block_modifiers.hh>
namespace ost{ namespace mol{ namespace mm{
std::vector<geom::Vec3> GromacsPositionRuleEvaluator::EvaluatePosRule(
int rule, int number,
const std::vector<geom::Vec3>& anchors) {
std::vector<geom::Vec3> positions;
switch(rule){
case 1:{
geom::Vec3 i_j = geom::Normalize(anchors[0]-anchors[1]);
geom::Vec3 i_k = geom::Normalize(anchors[0]-anchors[2]);
geom::Vec3 temp = geom::Normalize(i_j+i_k);
positions.push_back(anchors[0]+temp);
break;
}
case 2:{
Real x = 0.333806859234;
Real y = 0.942641491092;
geom::Vec3 i = anchors[0];
geom::Vec3 j = anchors[1];
geom::Vec3 k = anchors[2];
geom::Vec3 j_i = geom::Normalize(i-j);
geom::Vec3 q = geom::Normalize(geom::Cross(geom::Cross(k-j,j_i),j_i));
positions.push_back(i+x*j_i+y*q);
break;
}
case 3:{
Real x = 0.5;
Real y = 0.866025403784;
geom::Vec3 i = anchors[0];
geom::Vec3 j = anchors[1];
geom::Vec3 k = anchors[2];
geom::Vec3 j_i = geom::Normalize(i-j);
geom::Vec3 q = geom::Normalize(geom::Cross(geom::Cross(k-j,j_i),j_i));
positions.push_back(i+x*j_i-y*q);
positions.push_back(i+x*j_i+y*q);
break;
}
case 4:{
Real x = 0.333313247568;
Real y = 0.942816142732;
Real z = 0.866025403784;
geom::Vec3 i = anchors[0];
geom::Vec3 j = anchors[1];
geom::Vec3 k = anchors[2];
geom::Vec3 j_i = geom::Normalize(i-j);
geom::Vec3 q = geom::Normalize(geom::Cross(geom::Cross(k-j,j_i),j_i));
geom::Vec3 n1 = i+x*j_i+y*q;
geom::Vec3 temp = geom::Normalize(geom::Cross(j_i,n1-i));
geom::Vec3 n2 = i+x*j_i-0.5*q+z*temp;
geom::Vec3 n3 = i+x*j_i-0.5*q-z*temp;
positions.push_back(n1);
positions.push_back(n2);
if(number==3) positions.push_back(n3);
break;
}
case 5:{
geom::Vec3 i = anchors[0];
geom::Vec3 i_j = geom::Normalize(i-anchors[1]);
geom::Vec3 i_k = geom::Normalize(i-anchors[2]);
geom::Vec3 i_l = geom::Normalize(i-anchors[3]);
geom::Vec3 temp = geom::Normalize(i_j+i_k+i_l);
positions.push_back(i+temp);
break;
}
case 6:{
Real x = 0.81649043092;
Real y = 0.577358966516;
geom::Vec3 i = anchors[0];
geom::Vec3 i_j = geom::Normalize(i-anchors[1]);
geom::Vec3 i_k = geom::Normalize(i-anchors[2]);
geom::Vec3 temp_one = geom::Normalize(i_j+i_k);
geom::Vec3 temp_two = geom::Normalize(geom::Cross(i_j,i_k));
positions.push_back(i+y*temp_one+x*temp_two);
positions.push_back(i+y*temp_one-x*temp_two);
break;
}
case 7:{
//not implemented yet, do nothing
break;
}
case 8:{
Real x = 1.36*0.522498564716;
Real y = 1.36*0.852640164354;
geom::Vec3 i = anchors[0];
geom::Vec3 j = anchors[1];
geom::Vec3 k = anchors[2];
geom::Vec3 j_i = geom::Normalize(i-j);
geom::Vec3 q = geom::Normalize(geom::Cross(geom::Cross(k-j,j_i),j_i));
positions.push_back(i+x*j_i-y*q);
positions.push_back(i+x*j_i+y*q);
break;
}
case 9:{
Real x1 = 1.23*0.51503807491;
Real x2 = 1.25*0.422618261741;
Real y1 = 1.23*0.857167300702;
Real y2 = 1.25*0.906307787037;
geom::Vec3 i = anchors[0];
geom::Vec3 j = anchors[1];
geom::Vec3 k = anchors[2];
geom::Vec3 j_i = geom::Normalize(i-j);
geom::Vec3 q = geom::Normalize(geom::Cross(geom::Cross(k-j,j_i),j_i));
positions.push_back(i+x1*j_i-y1*q);
positions.push_back(i+x2*j_i+y2*q);
//we have to construct a third position hanging at the second one...
Real x = 0.333806859234;
Real y = 0.942641491092;
i = positions[1];
j = anchors[0];
k = anchors[1];
j_i = geom::Normalize(i-j);
q = geom::Normalize(geom::Cross(geom::Cross(k-j,j_i),j_i));
positions.push_back(i+x*j_i+y*q);
break;
}
case 10:{
//not implemented yet, do nothing
break;
}
case 11:{
//not implemented yet, do nothing
break;
}
default:{
break;
}
}
return positions;
}
void GromacsHydrogenConstructor::ApplyOnBuildingBlock(BuildingBlockPtr p) {
//we assume, that the hydrogens are already correct in the building block...
}
void GromacsHydrogenConstructor::ApplyOnResidue(ost::mol::ResidueHandle& res,
ost::mol::XCSEditor& ed) {
if(!res.IsValid()) return;
int number;
int method;
String prefix;
std::vector<String> anchor_names;
ost::mol::AtomHandleList anchor_atoms;
std::vector<String> hydrogen_names;
std::vector<geom::Vec3> hydrogen_positions;
std::vector<geom::Vec3> anchor_positions;
for(uint a=0;a<add_number_.size();++a){
hydrogen_names = hydrogen_names_[a];
number = add_number_[a];
method = methods_[a];
anchor_names = anchor_atom_names_[a];
hydrogen_positions.clear();
anchor_positions.clear();
anchor_atoms.clear();
String atom_name;
ost::mol::ResidueHandle temp_res;
ost::mol::AtomHandle atom;
bool found_anchors = true;
for (std::vector<String>::iterator it = anchor_names.begin();
it!=anchor_names.end(); ++it) {
temp_res = res;
atom_name = *it;
if(atom_name[0] == '-'){
atom_name.erase(0,1);
temp_res = res.GetPrev();
if(!temp_res.IsValid()){
//happens in case of termini... In this case the adding definitions
//are taken from the termini database...
found_anchors = false;
break;
}
}
if(atom_name[0] == '+'){
atom_name.erase(0,1);
temp_res = res.GetNext();
if(!temp_res.IsValid()){
//happens in case of termini... In this case the adding definitions
//are taken from the termini database...
found_anchors = false;
break;
}
}
atom = temp_res.FindAtom(atom_name);
if(!atom.IsValid()){
std::stringstream ss;
ss << "cannot find required anchor atom " << atom_name << " in residue ";
ss << res.GetQualifiedName() << " to add hydrogen!";
throw ost::Error(ss.str());
}
anchor_atoms.push_back(atom);
anchor_positions.push_back(atom.GetPos());
}
if(!found_anchors){
continue;
}
hydrogen_positions
= GromacsPositionRuleEvaluator::EvaluatePosRule(method, number,
anchor_positions);
for(int b=0;b<number;++b){
//only add hydrogen if not already present!
atom = res.FindAtom(hydrogen_names[b]);
if (!atom.IsValid()) {
atom = ed.InsertAtom(res, hydrogen_names[b], hydrogen_positions[b], "H");
}
else {
ed.SetAtomPos(atom, hydrogen_positions[b]);
}
}
}
}
void GromacsHydrogenConstructor::AddHydrogenRule(
uint number, int method,
const std::vector<String>& hydrogen_names,
const std::vector<String>& anchors) {
if(number != hydrogen_names.size()){
std::stringstream ss;
ss << "Tried to add hydrogen rule to hydrogen constructor with add number ";
ss << number << ", but "<< hydrogen_names.size() << "hydrogen names";
}
hydrogen_names_.push_back(hydrogen_names);
add_number_.push_back(number);
methods_.push_back(method);
anchor_atom_names_.push_back(anchors);
}
void GromacsBlockModifier::AddReplaceRule(const String& name, const String& new_name,
const String& new_type, Real new_charge){
replace_old_atom_name_.push_back(name);
replace_new_atom_name_.push_back(new_name);
replace_new_atom_type_.push_back(new_type);
replace_new_charge_.push_back(new_charge);
}
void GromacsBlockModifier::AddAddRule(int number, int method,
const std::vector<String>& atom_names,
const std::vector<String>& anchors,
const String& type, Real charge){
add_add_number_.push_back(number);
add_methods_.push_back(method);
add_atom_names_.push_back(atom_names);
add_anchor_atom_names_.push_back(anchors);
add_atom_types_.push_back(type);
add_charges_.push_back(charge);
}
void GromacsBlockModifier::ApplyOnBuildingBlock(BuildingBlockPtr p){
//let's resolve the replace statement
for(uint i = 0; i < replace_old_atom_name_.size(); ++i){
p->ReplaceAtom(replace_old_atom_name_[i], replace_new_atom_name_[i],
replace_new_atom_type_[i], replace_new_charge_[i]);
}
//let's resolve the add statement
for(uint i = 0; i < add_add_number_.size(); ++i){
String type = add_atom_types_[i];
Real charge = add_charges_[i];
std::vector<String> anchor_atoms = add_anchor_atom_names_[i];
std::vector<String> names = add_atom_names_[i];
//let's add the stuff to the buildingblock
for(std::vector<String>::iterator j = names.begin();
j != names.end(); ++j){
p->AddAtom(*j,type,charge);
//the bond, that connects the newly added atom to its anchor is not necessarily
//defined. Let's add it here. Note, that nothing happens, if this specific bond
//is already part of the building block.
InteractionPtr b_p(new Interaction(p->GetBonds()[0]->GetFuncType()));
std::vector<String> bond_names;
std::vector<String> bond_types;
bond_names.push_back(anchor_atoms[0]);
bond_names.push_back(*j);
bond_types.push_back(p->GetType(anchor_atoms[0]));
bond_types.push_back(type);
b_p->SetNames(bond_names);
b_p->SetTypes(bond_types);
p->AddBond(b_p);
}
}
//let's resolve the delete statement
for(std::vector<String>::iterator i = delete_atom_names_.begin();
i != delete_atom_names_.end(); ++i){
p->RemoveAtom(*i);
}
//we simply add the new interactions to the building block
for(std::vector<InteractionPtr>::iterator i = bonds_.begin();
i != bonds_.end(); ++i){
p->AddBond(*i, true);
}
for(std::vector<InteractionPtr>::iterator i = angles_.begin();
i != angles_.end(); ++i){
p->AddAngle(*i, true);
}
for(std::vector<InteractionPtr>::iterator i = dihedrals_.begin();
i != dihedrals_.end(); ++i){
p->AddDihedral(*i, true);
}
for(std::vector<InteractionPtr>::iterator i = impropers_.begin();
i != impropers_.end(); ++i){
p->AddImproper(*i, true);
}
for(std::vector<InteractionPtr>::iterator i = cmaps_.begin();
i != cmaps_.end(); ++i){
p->AddCMap(*i, true);
}
}
void GromacsBlockModifier::ApplyOnResidue(ost::mol::ResidueHandle& res,
ost::mol::XCSEditor& ed) {
//let's resolve the replace statement
for(uint i = 0; i < replace_old_atom_name_.size(); ++i){
ost::mol::AtomHandle at = res.FindAtom(replace_old_atom_name_[i]);
if(at.IsValid()){
ed.RenameAtom(at,replace_new_atom_name_[i]);
}
}
//let's resolve the add statement
for(uint i = 0; i < add_add_number_.size(); ++i){
int number = add_add_number_[i];
int method = add_methods_[i];
std::vector<String> anchor_names = add_anchor_atom_names_[i];
//String type = add_atom_types_[i];
std::vector<geom::Vec3> anchor_positions;
ost::mol::AtomHandleList anchor_atoms;
std::vector<String> names = add_atom_names_[i];
for(std::vector<String>::iterator j = anchor_names.begin();
j != anchor_names.end(); ++j){
ost::mol::AtomHandle at = res.FindAtom(*j);
if(!at.IsValid()){
std::stringstream ss;
ss << "Cannot find anchor atom required to construct termini";
ss << " in provided residue!";
throw ost::Error(ss.str());
}
anchor_positions.push_back(at.GetPos());
anchor_atoms.push_back(at);
}
std::vector<geom::Vec3> positions
= GromacsPositionRuleEvaluator::EvaluatePosRule(method, number,
anchor_positions);
String ele = "H";
if(method == 8 || method == 9){
ele = "O";
}
for(int j=0; j<number; ++j){
ost::mol::AtomHandle at = res.FindAtom(names[j]);
if(!at.IsValid()) at = ed.InsertAtom(res,names[j],positions[j],ele);
else ed.SetAtomPos(at,positions[j]);
ed.Connect(at, anchor_atoms[0]);
}
}
//let's resolve the delete statement
for(std::vector<String>::iterator i = delete_atom_names_.begin();
i != delete_atom_names_.end(); ++i){
ost::mol::AtomHandle at = res.FindAtom(*i);
if(at.IsValid()){
ed.DeleteAtom(at);
}
}
}
void GromacsBlockModifier::CheckInteractionToAdd(InteractionPtr p, const String& interaction_type) const{
if(p->GetNames().empty()) throw ost::Error("Expect interaction to have names properly set!");
if(interaction_type == "BOND"){
if(p->GetFuncType() != HarmonicBond){
throw ost::Error("Invalid function type encountered when adding bond!");
} return;
}
if(interaction_type == "ANGLE"){
if(p->GetFuncType() != HarmonicAngle && p->GetFuncType() != UreyBradleyAngle){
throw ost::Error("Invalid function type encountered when adding angle!");
} return;
}
if(interaction_type == "DIHEDRAL"){
if(p->GetFuncType() != PeriodicDihedral){
throw ost::Error("Invalid function type encountered when adding dihedral!");
} return;
}
if(interaction_type == "IMPROPER"){
if(p->GetFuncType() != PeriodicImproper && p->GetFuncType() != HarmonicImproper){
throw ost::Error("Invalid function type encountered when adding improper!");
} return;
}
if(interaction_type == "CMAP"){
if(p->GetFuncType() != CMap){
throw ost::Error("Invalid function type encountered when adding cmap!");
} return;
}
std::stringstream ss;
ss << "Cannot validate added interaction of type "<< interaction_type << "!";
throw ost::Error(ss.str());
}
}}}//ns