diff --git a/.gitignore b/.gitignore index 4378ad9..b302568 100644 --- a/.gitignore +++ b/.gitignore @@ -14,6 +14,14 @@ *.sty *.gif *.txt +*.crg +*.siz +*.cube +*.vscode +build/* +dist/* +*/talled +*.mp4 #Python *.pyc diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..097e9be --- /dev/null +++ b/.gitmodules @@ -0,0 +1,4 @@ +[submodule "cusplibrary"] + path = cusplibrary + url = https://github.com/lupoglaz/cusplibrary.git + branch = develop diff --git a/Benchmark/volumemul.py b/Benchmark/volumemul.py new file mode 100644 index 0000000..6fd3ec1 --- /dev/null +++ b/Benchmark/volumemul.py @@ -0,0 +1,64 @@ +from core import ModuleBenchmark +import torch +from torch.autograd import profiler +from TorchProteinLibrary import VolumeMul + +import numpy as np +import seaborn as sea +from matplotlib import pylab as plt +import pandas as pd +from tqdm import tqdm +import pickle as pkl + +class VolumeMulBenchmark(ModuleBenchmark): + def __init__(self, device='gpu', num_sequences=32, seq_length=350): + super().__init__(device) + self.vol1 = + + def prepare(self): + self.angles.requires_grad_() + self.coords = self.module(self.angles, self.length) + self.grad_coords = torch.ones(self.coords.size(), dtype=torch.float, device='cuda') + + def run_forward(self): + self.angles.detach() + self.coords = self.module(self.angles, self.length) + + def run_backward(self): + self.coords.backward(self.grad_coords) + +def test_length( interval=(60, 20, 300), num_measurements=1, output_filename='length.dat'): + time = [] + direction = [] + length = [] + for n in range(num_measurements): + for i in tqdm(range(interval[0], interval[2], interval[1])): + mb = Angles2BackboneBenchmark(seq_length=i) + dur_fwd = mb.measure_forward() + time.append(dur_fwd/1000.0) + direction.append('Forward') + length.append(i) + dur_bwd = mb.measure_backward() + time.append(dur_bwd/1000.0) + direction.append('Backward') + length.append(i) + + data = pd.DataFrame({ 'Time': time, + 'Direction': direction, + 'Length': length + }) + data.to_pickle(output_filename) + +if __name__=='__main__': + + test_length(interval=(60, 20, 700), output_filename='Data/ReducedModelTime.pkl', num_measurements=10) + + data = pd.read_pickle('Data/ReducedModelTime.pkl') + sea.set_style("whitegrid") + sea.set_context("paper", font_scale=1.5, rc={"lines.linewidth": 1.5}) + g1 = sea.relplot(x="Length", y="Time", hue='Direction', kind="line", style="Direction", height=6, aspect=1.3, markers=True, data=data) + plt.ylabel("Time, ms") + plt.xlabel("Sequence length, amino-acids") + sea.despine() + # plt.show() + plt.savefig("Fig/ReducedModelTime.png") \ No newline at end of file diff --git a/Layers/FullAtomModel/Angles2Coords/angles2coords_interface.cpp b/Layers/FullAtomModel/Angles2Coords/angles2coords_interface.cpp index b842ed6..5347547 100644 --- a/Layers/FullAtomModel/Angles2Coords/angles2coords_interface.cpp +++ b/Layers/FullAtomModel/Angles2Coords/angles2coords_interface.cpp @@ -9,7 +9,10 @@ void Angles2Coords_forward( torch::Tensor sequences, torch::Tensor input_angles, torch::Tensor output_coords, torch::Tensor res_names, - torch::Tensor atom_names + torch::Tensor res_nums, + torch::Tensor atom_names, + int polymer_type, + torch::Tensor chain_names ){ bool add_terminal = false; CHECK_CPU_INPUT_TYPE(sequences, torch::kByte); @@ -23,60 +26,122 @@ void Angles2Coords_forward( torch::Tensor sequences, } int batch_size = input_angles.sizes()[0]; - - #pragma omp parallel for - for(int i=0; i conf( seq, single_angles.data(), dummy_grad.data(), length, single_coords.data()); - for(int j=0; jatomNames.size(); k++){ - int idx = conf.groups[j]->atomIndexes[k]; - torch::Tensor single_atom_name = single_atom_names[idx]; - torch::Tensor single_res_name = single_res_names[idx]; - StringUtil::string2Tensor(ProtUtil::convertRes1to3(conf.groups[j]->residueName), single_res_name); - StringUtil::string2Tensor(conf.groups[j]->atomNames[k], single_atom_name); + + if(polymer_type == 0){ + + #pragma omp parallel for + + for(int i=0; i conf( seq, single_angles.data(), dummy_grad.data(), length, single_coords.data(), polymer_type, chain_names); + for(int j=0; jatomNames.size(); k++){ + int idx = conf.groups[j]->atomIndexes[k]; + torch::Tensor single_atom_name = single_atom_names[idx]; + torch::Tensor single_res_name = single_res_names[idx]; + single_res_nums[idx] = (int)conf.groups[j]->residueIndex; + StringUtil::string2Tensor(ProtUtil::convertRes1to3(conf.groups[j]->residueName), single_res_name); + StringUtil::string2Tensor(conf.groups[j]->atomNames[k], single_atom_name); + } } - } - })); - // cConformation conf( seq, single_angles.data(), dummy_grad.data(), - // length, single_coords.data()); - //Output atom names and residue names - + })); + // cConformation conf( seq, single_angles.data(), dummy_grad.data(), + // length, single_coords.data()); + //Output atom names and residue names + } + } + if(polymer_type == 1 or polymer_type == 2){ + + #pragma omp parallel for + + for(int i=0; i conf( seq, single_angles.data(), dummy_grad.data(), length, single_coords.data(), polymer_type, chain_names); + for(int j=0; jatomNames.size(); k++){ + int idx = conf.groups[j]->atomIndexes[k]; + torch::Tensor single_atom_name = single_atom_names[idx]; + torch::Tensor single_res_name = single_res_names[idx]; + single_res_nums[idx] = (int)conf.groups[j]->residueIndex; + StringUtil::string2Tensor(ProtUtil::convertRes1to3(conf.groups[j]->residueName, polymer_type), single_res_name); + StringUtil::string2Tensor(conf.groups[j]->atomNames[k], single_atom_name);} + } + })); + } } } void Angles2Coords_backward( torch::Tensor grad_atoms, torch::Tensor grad_angles, torch::Tensor sequences, - torch::Tensor input_angles + torch::Tensor input_angles, + int polymer_type, + torch::Tensor chain_names ){ bool add_terminal = false; CHECK_CPU_INPUT_TYPE(sequences, torch::kByte); @@ -98,13 +163,15 @@ void Angles2Coords_backward( torch::Tensor grad_atoms, torch::Tensor single_grad_atoms = grad_atoms[i]; std::string seq = StringUtil::tensor2String(single_sequence); - + +// Get Num Atoms uint length = single_angles.sizes()[1]; - int num_atoms = ProtUtil::getNumAtoms(seq, add_terminal); + int num_atoms = ProtUtil::getNumAtoms(seq, add_terminal, polymer_type, chain_names); +// Conformation torch::Tensor dummy_coords = torch::zeros({3*num_atoms}, torch::TensorOptions().dtype(grad_atoms.dtype())); AT_DISPATCH_FLOATING_TYPES(single_angles.type(), "cConformation", ([&] { - cConformation conf(seq, single_angles.data(), single_grad_angles.data(),length, dummy_coords.data()); + cConformation conf(seq, single_angles.data(), single_grad_angles.data(),length, dummy_coords.data(), polymer_type, chain_names); conf.backward(conf.root, single_grad_atoms.data()); })); // cConformation conf( seq, single_angles.data(), single_grad_angles.data(), @@ -125,7 +192,8 @@ void Angles2Coords_save( const char* sequence, } std::string aa(sequence); uint length = aa.length(); - int num_atoms = ProtUtil::getNumAtoms(aa, add_terminal); + int polymer_type = 1; //"hard coded" need to change to arg + int num_atoms = ProtUtil::getNumAtoms(aa, add_terminal, polymer_type); torch::Tensor dummy_grad = torch::zeros_like(input_angles); torch::Tensor dummy_coords = torch::zeros({3*num_atoms}, torch::TensorOptions().dtype(input_angles.dtype())); AT_DISPATCH_FLOATING_TYPES(input_angles.type(), "cConformation.save", ([&]{ @@ -136,9 +204,9 @@ void Angles2Coords_save( const char* sequence, } -int getSeqNumAtoms( const char *sequence){ +int getSeqNumAtoms( const char *sequence, int polymer_type, torch::Tensor chain_names){ bool add_terminal = false; std::string seq(sequence); - int num_atoms = ProtUtil::getNumAtoms(seq, add_terminal); + int num_atoms = ProtUtil::getNumAtoms(seq, add_terminal, polymer_type, chain_names); return num_atoms; } diff --git a/Layers/FullAtomModel/Angles2Coords/angles2coords_interface.h b/Layers/FullAtomModel/Angles2Coords/angles2coords_interface.h index 3c1a79a..fd1deed 100644 --- a/Layers/FullAtomModel/Angles2Coords/angles2coords_interface.h +++ b/Layers/FullAtomModel/Angles2Coords/angles2coords_interface.h @@ -3,17 +3,22 @@ void Angles2Coords_forward( at::Tensor sequences, at::Tensor input_angles, at::Tensor output_coords, at::Tensor res_names, - at::Tensor atom_names + torch::Tensor res_nums, + at::Tensor atom_names, + int polymer_type, + torch::Tensor chain_names ); void Angles2Coords_backward( at::Tensor grad_atoms, at::Tensor grad_angles, at::Tensor sequences, - at::Tensor input_angles + at::Tensor input_angles, + int polymer_type, + torch::Tensor chain_names ); void Angles2Coords_save( const char* sequence, at::Tensor input_angles, const char* output_filename, const char mode ); -int getSeqNumAtoms( const char *sequence); \ No newline at end of file +int getSeqNumAtoms( const char *sequence, int polymer_type = 0, torch::Tensor chain_names = {}); \ No newline at end of file diff --git a/Layers/FullAtomModel/Coords2TypedCoords/coords2typedcoords_interface.cpp b/Layers/FullAtomModel/Coords2TypedCoords/coords2typedcoords_interface.cpp index 483f78e..74cf847 100644 --- a/Layers/FullAtomModel/Coords2TypedCoords/coords2typedcoords_interface.cpp +++ b/Layers/FullAtomModel/Coords2TypedCoords/coords2typedcoords_interface.cpp @@ -11,20 +11,19 @@ void Coords2TypedCoords_forward( torch::Tensor input_coords, torch::Tensor input_num_atoms, torch::Tensor output_coords, torch::Tensor output_num_atoms_of_type, - torch::Tensor output_offsets, - torch::Tensor output_atom_indexes + torch::Tensor output_atom_indexes, + int num_atom_types ){ - const uint num_atom_types = 11; + //const uint num_atom_types = 4; CHECK_CPU_INPUT(input_coords); CHECK_CPU_INPUT_TYPE(res_names, torch::kByte); CHECK_CPU_INPUT_TYPE(atom_names, torch::kByte); CHECK_CPU_INPUT_TYPE(input_num_atoms, torch::kInt); CHECK_CPU_INPUT(output_coords); CHECK_CPU_INPUT_TYPE(output_num_atoms_of_type, torch::kInt); - CHECK_CPU_INPUT_TYPE(output_offsets, torch::kInt); CHECK_CPU_INPUT_TYPE(output_atom_indexes, torch::kInt); - if(input_coords.ndimension() != 2){ + if(input_coords.ndimension() != 3){ ERROR("Incorrect input ndim"); } int batch_size = input_coords.size(0); @@ -33,58 +32,44 @@ void Coords2TypedCoords_forward( torch::Tensor input_coords, for(int i=0; i()[i]; - torch::Tensor single_intput_coords = input_coords[i]; + torch::Tensor single_input_coords = input_coords[i]; torch::Tensor single_atom_names = atom_names[i]; torch::Tensor single_res_names = res_names[i]; torch::Tensor single_output_coords = output_coords[i]; torch::Tensor single_output_num_atoms_of_type = output_num_atoms_of_type[i]; - torch::Tensor single_output_offsets = output_offsets[i]; torch::Tensor single_output_atom_indexes = output_atom_indexes[i]; - int num_atoms_added[num_atom_types]; - int atom_types[num_atoms]; - for(int j=0;j(); + auto a_num_atoms_of_type = single_output_num_atoms_of_type.accessor(); - //Assign atom types - torch::Tensor single_atom_name, single_res_name; - + int type; for(int j=0; j(); - for(int j=0;j r_dst( single_output_coords.data() + 3*dst_idx ); - cVector3 r_src( single_intput_coords.data() + 3*j); + cVector3 r_dst( single_output_coords[type].data() + 3*dst_idx ); + cVector3 r_src( single_input_coords.data() + 3*j); r_dst = r_src; })); - single_output_atom_indexes[offset + num_atoms_added[type]] = j; - num_atoms_added[type] += 1; + a_atom_indexes[type][dst_idx] = j; + a_num_atoms_of_type[type] += 1; } } @@ -93,17 +78,16 @@ void Coords2TypedCoords_forward( torch::Tensor input_coords, void Coords2TypedCoords_backward( torch::Tensor grad_typed_coords, torch::Tensor grad_flat_coords, torch::Tensor num_atoms_of_type, - torch::Tensor offsets, - torch::Tensor atom_indexes + torch::Tensor atom_indexes, + int num_atom_types ){ - const uint num_atom_types=11; + //const uint num_atom_types=4; CHECK_CPU_INPUT(grad_typed_coords); CHECK_CPU_INPUT(grad_flat_coords); CHECK_CPU_INPUT_TYPE(num_atoms_of_type, torch::kInt); - CHECK_CPU_INPUT_TYPE(offsets, torch::kInt); CHECK_CPU_INPUT_TYPE(atom_indexes, torch::kInt); - if(grad_typed_coords.ndimension() != 2){ + if(grad_typed_coords.ndimension() != 3){ ERROR("Incorrect input ndim"); } @@ -114,21 +98,20 @@ void Coords2TypedCoords_backward( torch::Tensor grad_typed_coords, torch::Tensor single_grad_typed_coords = grad_typed_coords[i]; torch::Tensor single_grad_flat_coords = grad_flat_coords[i]; torch::Tensor single_atom_indexes = atom_indexes[i]; - torch::Tensor single_offsets = offsets[i]; torch::Tensor single_num_atoms_of_type = num_atoms_of_type[i]; - for(int j=0;j()[j]; - int offset = single_offsets.accessor()[j]; - for(int k=0; k()[src_idx]; + auto a_single_atom_indexes = single_atom_indexes.accessor(); + auto a_single_num_atoms_of_type = single_num_atoms_of_type.accessor(); + + for(int type=0; type r_src( single_grad_typed_coords.data()+ 3*src_idx ); - cVector3 r_dst( single_grad_flat_coords.data() + 3*dst_idx ); + cVector3 r_src( single_grad_typed_coords[type].data()+ 3*at_idx ); + cVector3 r_dst( single_grad_flat_coords.data() + 3*a_single_atom_indexes[type][at_idx] ); r_dst = r_src; })); } } } -} \ No newline at end of file +} diff --git a/Layers/FullAtomModel/Coords2TypedCoords/coords2typedcoords_interface.h b/Layers/FullAtomModel/Coords2TypedCoords/coords2typedcoords_interface.h index 55af9b1..d00b886 100644 --- a/Layers/FullAtomModel/Coords2TypedCoords/coords2typedcoords_interface.h +++ b/Layers/FullAtomModel/Coords2TypedCoords/coords2typedcoords_interface.h @@ -5,12 +5,12 @@ void Coords2TypedCoords_forward( torch::Tensor input_coords, torch::Tensor input_num_atoms, torch::Tensor output_coords, torch::Tensor output_num_atoms_of_type, - torch::Tensor output_offsets, - torch::Tensor output_atom_indexes + torch::Tensor output_atom_indexes, + int num_atom_types ); void Coords2TypedCoords_backward( torch::Tensor grad_typed_coords, torch::Tensor grad_flat_coords, torch::Tensor num_atoms_of_type, - torch::Tensor offsets, - torch::Tensor atom_indexes - ); \ No newline at end of file + torch::Tensor atom_indexes, + int num_atom_types + ); diff --git a/Layers/FullAtomModel/CoordsTransform/coordsTransform_interface.cpp b/Layers/FullAtomModel/CoordsTransform/coordsTransform_interface.cpp index 2497018..1a33e10 100644 --- a/Layers/FullAtomModel/CoordsTransform/coordsTransform_interface.cpp +++ b/Layers/FullAtomModel/CoordsTransform/coordsTransform_interface.cpp @@ -83,6 +83,7 @@ void CoordsRotate_forward( torch::Tensor input_coords, AT_DISPATCH_FLOATING_TYPES(input_coords.type(), "CoordsRotate_forward", ([&]{ cMatrix33 _R = tensor2Matrix33(single_R); + rotate(single_input_coords, _R, single_output_coords, num_at[i]); })); } diff --git a/Layers/FullAtomModel/PDB2Coords/pdb2coords_interface.cpp b/Layers/FullAtomModel/PDB2Coords/pdb2coords_interface.cpp index 297f7d8..345a18a 100644 --- a/Layers/FullAtomModel/PDB2Coords/pdb2coords_interface.cpp +++ b/Layers/FullAtomModel/PDB2Coords/pdb2coords_interface.cpp @@ -4,69 +4,245 @@ #include #include #include -/* -void PDB2CoordsOrdered(torch::Tensor filenames, torch::Tensor coords, torch::Tensor res_names, torch::Tensor atom_names, torch::Tensor num_atoms, torch::Tensor mask){ - bool add_terminal = true; - if( filenames.dtype() != torch::kByte || res_names.dtype() != torch::kByte || atom_names.dtype() != torch::kByte - || coords.dtype() != torch::kDouble || num_atoms.dtype() != torch::kInt || mask.dtype() != torch::kByte){ - throw("Incorrect tensor types"); - } - if(filenames.ndimension() != 2){ - throw("Incorrect input ndim"); + +void PDB2CoordsOrdered( torch::Tensor filenames, torch::Tensor coords, torch::Tensor chain_names, torch::Tensor res_names, + torch::Tensor res_nums, torch::Tensor atom_names, torch::Tensor atom_mask, torch::Tensor num_atoms, torch::Tensor chain_ids, int polymer_type){ + CHECK_CPU_INPUT_TYPE(filenames, torch::kByte); + CHECK_CPU_INPUT_TYPE(res_names, torch::kByte); + CHECK_CPU_INPUT_TYPE(atom_names, torch::kByte); + CHECK_CPU_INPUT_TYPE(chain_names, torch::kByte); + CHECK_CPU_INPUT_TYPE(res_nums, torch::kInt); + CHECK_CPU_INPUT_TYPE(num_atoms, torch::kInt); + CHECK_CPU_INPUT_TYPE(atom_mask, torch::kByte); + CHECK_CPU_INPUT(coords); + if(coords.ndimension() != 2){ + ERROR("Incorrect input ndim"); } + if (polymer_type == 0){ + int batch_size = filenames.size(0); + std::string resLastAtom("OXT"); + #pragma omp parallel for + for(int i=0; i 0){ + torch::Tensor single_chain_id = chain_ids[i]; + chain_id = StringUtil::tensor2String(single_chain_id); + } + cPDBLoader pdb(filename, chain_id, 0); + num_atoms[i] = 0; + int previous_res_num = pdb.res_nums[0]; + for(int j=0; j()[0]; + int64_t size_nums[] = {batch_size, max_num_atoms}; + int64_t size_coords[] = {batch_size, max_num_atoms*3}; + int64_t size_names[] = {batch_size, max_num_atoms, 4}; + + coords.resize_(torch::IntList(size_coords, 2)).fill_(0.0); + chain_names.resize_(torch::IntList(size_names, 3)).fill_(0); + res_names.resize_(torch::IntList(size_names, 3)).fill_(0); + res_nums.resize_(torch::IntList(size_nums, 2)).fill_(0); + atom_names.resize_(torch::IntList(size_names, 3)).fill_(0); + atom_mask.resize_(torch::IntList(size_nums, 2)).fill_(0); + + + #pragma omp parallel for + for(int i=0; i 0){ + torch::Tensor single_chain_id = chain_ids[i]; + chain_id = StringUtil::tensor2String(single_chain_id); + } + + cPDBLoader pdb(filename, chain_id, 0); + + int global_ind = 0; + int previous_res_num = pdb.res_nums[0]; + for(int j=0; j()[0]; - - int64_t size_coords[] = {batch_size, max_num_atoms*3}; - int64_t size_names[] = {batch_size, max_num_atoms, 4}; - int64_t size_mask[] = {batch_size, max_num_atoms}; - - coords.resize_(torch::IntList(size_coords, 2)); - res_names.resize_(torch::IntList(size_names, 3)); - atom_names.resize_(torch::IntList(size_names, 3)); - mask.resize_(torch::IntList(size_mask, 3)); - - #pragma omp parallel for - for(int i=0; i 0){ + torch::Tensor single_chain_id = chain_ids[i]; + chain_id = StringUtil::tensor2String(single_chain_id); + + cPDBLoader pdb(filename, chain_id, polymer_type); + + num_atoms[i] = 0; + + + num_atoms[i] = static_cast(pdb.atom_names.size() + 1); + +// Get Num_atoms another way : I can use the following way now, I modified getAtomIndex just need to add the correct args +// Then I might be able to combine this first part of poly 0 & 1+2. (could i just use getNumAtom?) +// Then once I get poly type I can use an if to control flow of parts of the 2nd part I need to + +// int previous_res_num = pdb.res_nums[0]; +// for(int j=0; j()[0]; +// int max_num_atoms = num_atoms + int64_t size_nums[] = {batch_size, max_num_atoms}; + int64_t size_coords[] = {batch_size, max_num_atoms*3}; + int64_t size_names[] = {batch_size, max_num_atoms, 4}; +// + coords.resize_(torch::IntList(size_coords, 2)).fill_(0.0); + chain_names.resize_(torch::IntList(size_names, 3)).fill_(0); + res_names.resize_(torch::IntList(size_names, 3)).fill_(0); + res_nums.resize_(torch::IntList(size_nums, 2)).fill_(0); + atom_names.resize_(torch::IntList(size_names, 3)).fill_(0); + atom_mask.resize_(torch::IntList(size_nums, 2)).fill_(0); +// +// + #pragma omp parallel for + for(int i=0; i 0){ + torch::Tensor single_chain_id = chain_ids[i]; + chain_id = StringUtil::tensor2String(single_chain_id); } - global_ind += ProtUtil::getAtomIndex(pdb.res_res_names[j], lastO) + 1; + cPDBLoader pdb(filename, chain_id, polymer_type); + + int global_ind = 0; + int previous_res_num = pdb.res_nums[0]; + std::string chain_idx; + int chain_num = 0; + + for(int j=0; j chain_idx && pdb.atom_names[j] == "O5'"){ + chain_idx = pdb.chain_names[j]; + int res_idx = static_cast(pdb.res_nums[j]); + + while (pdb.res_nums[j] == pdb.res_nums[res_idx]){ + bool fiveprime_ind = true; + + if (previous_res_num < pdb.res_nums[j]) { + previous_res_num = pdb.res_nums[j]; + if (pdb.res_names[j-1] == "DA" || pdb.res_names[j-1] == "DG" || pdb.res_names[j-1] == "A" || pdb.res_names[j-1] == "G") { + std::string resLastAtom("C4"); + global_ind += ProtUtil::getAtomIndex(pdb.res_names[j-1], resLastAtom, true, polymer_type) + 1; + } + if (pdb.res_names[j-1] == "DT" || pdb.res_names[j-1] == "DC" || pdb.res_names[j-1] == "C" || pdb.res_names[j-1] == "U") { + std::string resLastAtom("C6"); + global_ind += ProtUtil::getAtomIndex(pdb.res_names[j-1], resLastAtom, true, polymer_type) + 1; + } + } + uint idx = ProtUtil::getAtomIndex(pdb.res_names[j], pdb.atom_names[j], true, polymer_type) + global_ind - (3 * chain_num); + + StringUtil::string2Tensor(pdb.chain_names[j], single_chain_names[idx]); + StringUtil::string2Tensor(pdb.res_names[j], single_res_names[idx]); + StringUtil::string2Tensor(pdb.atom_names[j], single_atom_names[idx]); + single_res_nums[idx] = pdb.res_nums[j]; + + single_coords[3*idx + 0] = pdb.r[j].v[0]; + single_coords[3*idx + 1] = pdb.r[j].v[1]; + single_coords[3*idx + 2] = pdb.r[j].v[2]; + single_mask[idx] = 1; + + ++j; + } + ++chain_num; + } + if (previous_res_num < pdb.res_nums[j]) { + previous_res_num = pdb.res_nums[j]; + if (pdb.res_names[j-1] == "DA" || pdb.res_names[j-1] == "DG" || pdb.res_names[j-1] == "A" || pdb.res_names[j-1] == "G") { + std::string resLastAtom("C4"); + global_ind += ProtUtil::getAtomIndex(pdb.res_names[j-1], resLastAtom, false, polymer_type) + 1; + } + if (pdb.res_names[j-1] == "DT" || pdb.res_names[j-1] == "DC" || pdb.res_names[j-1] == "C" || pdb.res_names[j-1] == "U") { + std::string resLastAtom("C6"); + global_ind += ProtUtil::getAtomIndex(pdb.res_names[j-1], resLastAtom, false, polymer_type) + 1; + } + } + uint idx = ProtUtil::getAtomIndex(pdb.res_names[j], pdb.atom_names[j], false, polymer_type) + global_ind - (3 * chain_num); + + + StringUtil::string2Tensor(pdb.chain_names[j], single_chain_names[idx]); + StringUtil::string2Tensor(pdb.res_names[j], single_res_names[idx]); + StringUtil::string2Tensor(pdb.atom_names[j], single_atom_names[idx]); + single_res_nums[idx] = pdb.res_nums[j]; + + single_coords[3*idx + 0] = pdb.r[j].v[0]; + single_coords[3*idx + 1] = pdb.r[j].v[1]; + single_coords[3*idx + 2] = pdb.r[j].v[2]; + single_mask[idx] = 1; + } +// } +// } } } + + else{ + std::cerr << "Error Polymer Type Is Not Valid \n"; + } } -*/ + void PDB2CoordsUnordered( torch::Tensor filenames, torch::Tensor coords, torch::Tensor chain_names, torch::Tensor res_names, torch::Tensor res_nums, torch::Tensor atom_names, torch::Tensor num_atoms){ CHECK_CPU_INPUT_TYPE(filenames, torch::kByte); @@ -86,7 +262,7 @@ void PDB2CoordsUnordered( torch::Tensor filenames, torch::Tensor coords, torch for(int i=0; i()[0]; @@ -94,11 +270,11 @@ void PDB2CoordsUnordered( torch::Tensor filenames, torch::Tensor coords, torch int64_t size_coords[] = {batch_size, max_num_atoms*3}; int64_t size_names[] = {batch_size, max_num_atoms, 4}; - coords.resize_(torch::IntList(size_coords, 2)); - chain_names.resize_(torch::IntList(size_names, 3)); - res_names.resize_(torch::IntList(size_names, 3)); - res_nums.resize_(torch::IntList(size_nums, 2)); - atom_names.resize_(torch::IntList(size_names, 3)); + coords.resize_(torch::IntList(size_coords, 2)).fill_(0); + chain_names.resize_(torch::IntList(size_names, 3)).fill_(0); + res_names.resize_(torch::IntList(size_names, 3)).fill_(0); + res_nums.resize_(torch::IntList(size_nums, 2)).fill_(0); + atom_names.resize_(torch::IntList(size_names, 3)).fill_(0); #pragma omp parallel for for(int i=0; i r_target(single_coords.data() + 3*j); diff --git a/Layers/FullAtomModel/PDB2Coords/pdb2coords_interface.h b/Layers/FullAtomModel/PDB2Coords/pdb2coords_interface.h index ebbc78c..aebf709 100644 --- a/Layers/FullAtomModel/PDB2Coords/pdb2coords_interface.h +++ b/Layers/FullAtomModel/PDB2Coords/pdb2coords_interface.h @@ -1,3 +1,4 @@ #include -// void PDB2CoordsOrdered(torch::Tensor filenames, torch::Tensor coords, torch::Tensor res_names, torch::Tensor atom_names, torch::Tensor mask); +void PDB2CoordsOrdered( torch::Tensor filenames, torch::Tensor coords, torch::Tensor chain_names, torch::Tensor res_names, + torch::Tensor res_nums, torch::Tensor atom_names, torch::Tensor atom_mask, torch::Tensor num_atoms, torch::Tensor chain_ids, int polymer_type); void PDB2CoordsUnordered(torch::Tensor filenames, torch::Tensor coords, torch::Tensor chain_names, torch::Tensor res_names, torch::Tensor res_nums, torch::Tensor atom_names, torch::Tensor num_atoms); diff --git a/Layers/FullAtomModel/cConformation.cpp b/Layers/FullAtomModel/cConformation.cpp index 75cb212..2d90a0a 100644 --- a/Layers/FullAtomModel/cConformation.cpp +++ b/Layers/FullAtomModel/cConformation.cpp @@ -9,7 +9,7 @@ template void cTransform::updateMatrix(){ Ry.setRy(*beta); Tr.setT(d, 'x'); Rx.setRx(*alpha); - mat = Ry*Tr*Rx; + mat = Ry*Tr*Rx; // Rx is dihedral angles (rotation/transform this) } template void cTransform::updateDMatrix(){ cMatrix44 Ry, Tr, DRx; @@ -26,12 +26,12 @@ template void cTransform::print(){ // template std::ostream& operator<<(std::ostream& os, const cNode& node){ // return os<<*(node.group); // } - -template cConformation::cConformation(std::string aa, T *angles, T *angles_grad, uint angles_length, T *atoms_global, bool add_terminal){ +template cConformation::cConformation(std::string aa, T *angles, T *angles_grad, uint angles_length, T *atoms_global, int polymer_type, torch::Tensor chain_names,bool add_terminal){ cNode *lastC = NULL; zero_const = 0.0; this->atoms_global = atoms_global; bool terminal = false; + if( polymer_type == 0){ for(int i=0; i cConformation::cConformation(std::string aa, T *angles, break; } } - + } + if( polymer_type == 1){ + std::string chain_idx = "0"; + torch::Tensor single_chain_names = chain_names[0]; + int atom_idx = 0; + char last_res = '0'; + + for(int i=0; i params({alpha, beta, gamma, delta, epsilon, zeta, nu0, nu1, nu2, nu3, nu4, chi}); + std::vector params_grad({dalpha, dbeta, dgamma, ddelta, depsilon, dzeta, dnu0, dnu1, dnu2, dnu3, dnu4, dchi}); + + + if(tensor2String(single_chain_names[atom_idx]) > chain_idx){ + chain_idx = tensor2String(single_chain_names[atom_idx]); + if (aa[i] == 'G'){ + lastC = addDG_5Prime(lastC, params, params_grad, last_res); + std::string term_atom("C4"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, true, polymer_type) + 1; + last_res = 'G'; + continue; + } + if (aa[i] == 'A'){ + lastC = addDA_5Prime(lastC, params, params_grad, last_res); + std::string term_atom("C4"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, true, polymer_type) + 1; + last_res = 'A'; + continue; + } + if (aa[i] == 'T'){ + lastC = addDT_5Prime(lastC, params, params_grad, last_res); + std::string term_atom("C6"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, true, polymer_type) + 1; + last_res = 'T'; + continue; + } + if (aa[i] == 'C'){ + lastC = addDC_5Prime(lastC, params, params_grad, last_res); + std::string term_atom("C6"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, true, polymer_type) + 1; + last_res = 'C'; + continue; + } + } + + if (aa[i] == 'G'){ + lastC = addDG(lastC, params, params_grad, last_res); + std::string term_atom("C4"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, false, polymer_type) + 1; + last_res = 'G'; + } + if (aa[i] == 'A'){ + lastC = addDA(lastC, params, params_grad, last_res); //addDA + std::string term_atom("C4"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, false, polymer_type) + 1; + last_res = 'A'; + } + if (aa[i] == 'T'){ + lastC = addDT(lastC, params, params_grad, last_res); //addDT + std::string term_atom("C6"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, false, polymer_type) + 1; + last_res = 'T'; + } + if (aa[i] == 'C'){ + lastC = addDC(lastC, params, params_grad, last_res); //addDC + std::string term_atom("C6"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, false, polymer_type) + 1; + last_res = 'C'; + } + } + } + if (polymer_type == 2){ + std::string chain_idx = "0"; + torch::Tensor single_chain_names = chain_names[0]; + int atom_idx = 0; + char last_res = '0'; + + for(int i=0; i params({alpha, beta, gamma, delta, epsilon, zeta, nu0, nu1, nu2, nu3, nu4, chi}); + std::vector params_grad({dalpha, dbeta, dgamma, ddelta, depsilon, dzeta, dnu0, dnu1, dnu2, dnu3, dnu4, dchi}); + + if(tensor2String(single_chain_names[atom_idx]) > chain_idx){ + chain_idx = tensor2String(single_chain_names[atom_idx]); + if (aa[i] == 'G'){ + lastC = addG_5Prime(lastC, params, params_grad, last_res); + std::string term_atom("C4"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, true, polymer_type) + 1; + last_res = 'G'; + continue; + } + if (aa[i] == 'A'){ + lastC = addA_5Prime(lastC, params, params_grad, last_res); //addDA + std::string term_atom("C4"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, true, polymer_type) + 1; + last_res = 'A'; + continue; + } + if (aa[i] == 'U'){ + lastC = addU_5Prime(lastC, params, params_grad, last_res); //addDT + std::string term_atom("C6"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, true, polymer_type) + 1; + last_res = 'U'; + continue; + } + if (aa[i] == 'C'){ + lastC = addC_5Prime(lastC, params, params_grad, last_res); //addDC + std::string term_atom("C6"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, true, polymer_type) + 1; + last_res = 'C'; + continue; + } + } + + if (aa[i] == 'G'){ + lastC = addG(lastC, params, params_grad, last_res); //addDG *terminal == five_prime + std::string term_atom("C4"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, false, polymer_type) + 1; + last_res = 'G'; + } + if (aa[i] == 'A'){ + lastC = addA(lastC, params, params_grad, last_res); //addDA + std::string term_atom("C4"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, false, polymer_type) + 1; + last_res = 'A'; + } + if (aa[i] == 'U'){ + lastC = addU(lastC, params, params_grad, last_res); //addDT + std::string term_atom("C6"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, false, polymer_type) + 1; + last_res = 'U'; + } + if (aa[i] == 'C'){ + lastC = addC(lastC, params, params_grad, last_res); //addDC + std::string term_atom("C6"); + std::string NA(1, aa[i]); + atom_idx += getAtomIndex(NA, term_atom, false, polymer_type) + 1; + last_res = 'C'; + } + } + } + + //Computing conformation this->update(this->root); //Computing number of atoms @@ -138,7 +319,6 @@ template cConformation::~cConformation(){ for(int i=0; i cNode *cConformation::addNode(cNode *parent, cRigidGroup *group, cTransform *t){ cNode *new_node = new cNode(); new_node->group = group; @@ -224,7 +404,6 @@ template void cConformation::print(cNode *node){ // std::cout<<".\n"; } - template void cConformation::save(std::string filename, const char mode){ if(mode=='w'){ std::ofstream pfile(filename, std::ofstream::out); diff --git a/Layers/FullAtomModel/cConformation.h b/Layers/FullAtomModel/cConformation.h index 08405da..38c1701 100644 --- a/Layers/FullAtomModel/cConformation.h +++ b/Layers/FullAtomModel/cConformation.h @@ -6,7 +6,8 @@ #include #include #include - +#include +//AS: function for graphs template class cTransform{ public: @@ -62,7 +63,7 @@ class cConformation{ cNode *root; // Construct protein graph and bind grad to the angles - cConformation(std::string aa, T *angles, T *angles_grad, uint angles_length, T *atoms_global, bool add_terminal=false); + cConformation(std::string aa, T *angles, T *angles_grad, uint angles_length, T *atoms_global, int polymer_type = 0, torch::Tensor chain_names = {}, bool add_terminal=false); ~cConformation(); void print(cNode *node); @@ -100,6 +101,22 @@ class cConformation{ cNode *addPhe(cNode *parentC, std::vector params, std::vector params_grad, bool terminal=false); cNode *addTyr(cNode *parentC, std::vector params, std::vector params_grad, bool terminal=false); cNode *addTrp(cNode *parentC, std::vector params, std::vector params_grad, bool terminal=false); + cNode *addDG(cNode *parentC, std::vector params, std::vector params_grad, char last_res); + cNode *addDA(cNode *parentC, std::vector params, std::vector params_grad, char last_res); + cNode *addDT(cNode *parentC, std::vector params, std::vector params_grad, char last_res); + cNode *addDC(cNode *parentC, std::vector params, std::vector params_grad, char last_res); + cNode *addDG_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res); + cNode *addDA_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res); + cNode *addDT_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res); + cNode *addDC_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res); + cNode *addG(cNode *parentC, std::vector params, std::vector params_grad, char last_res); + cNode *addA(cNode *parentC, std::vector params, std::vector params_grad, char last_res); + cNode *addU(cNode *parentC, std::vector params, std::vector params_grad, char last_res); + cNode *addC(cNode *parentC, std::vector params, std::vector params_grad, char last_res); + cNode *addG_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res); + cNode *addA_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res); + cNode *addU_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res); + cNode *addC_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res); }; #endif \ No newline at end of file diff --git a/Layers/FullAtomModel/cConformationAA.cpp b/Layers/FullAtomModel/cConformationAA.cpp index 9e282ab..7b565e2 100644 --- a/Layers/FullAtomModel/cConformationAA.cpp +++ b/Layers/FullAtomModel/cConformationAA.cpp @@ -1,4 +1,5 @@ #include "cConformation.h" +#include #define PARENT_CHECK \ if(parentC!=NULL){ \ @@ -7,7 +8,53 @@ }else{ \ residueIndex = 0; \ firstAtomIndex = 0; \ - } + } + +#define PARENT_CHECK_DNA \ + if(parentC!=NULL){ \ + int shift = 8;\ + switch(last_res){\ + case 'C': \ + shift = 11;\ + break;\ + case 'T': \ + shift = 12;\ + break;\ + case 'A': \ + shift = 13;\ + break;\ + case 'G': \ + shift = 14;\ + break;\ + }residueIndex = parentC->group->residueIndex + 1; \ + firstAtomIndex = parentC->group->atomIndexes.back() + shift;\ + }else{ \ + residueIndex = 0; \ + firstAtomIndex = 0; \ + } + +#define PARENT_CHECK_RNA \ + if(parentC!=NULL){ \ + int shift = 8;\ + switch(last_res){\ + case 'C': \ + shift = 12;\ + break;\ + case 'U': \ + shift = 12;\ + break;\ + case 'A': \ + shift = 14;\ + break;\ + case 'G': \ + shift = 15;\ + break;\ + }residueIndex = parentC->group->residueIndex + 1; \ + firstAtomIndex = parentC->group->atomIndexes.back() + shift;\ + }else{ \ + residueIndex = 0; \ + firstAtomIndex = 0; \ + } #define ADD_NITROGEN \ bbN = makeAtom("N", firstAtomIndex, residueName, residueIndex, atoms_global); \ @@ -61,6 +108,170 @@ this->transforms.push_back(bbC_transform); \ nC = addNode(nCA, groups.back(), transforms.back()); +#define ADD_PHOSPHATE \ + bbP = makePhosGroup(geo, firstAtomIndex, residueName, residueIndex, atoms_global); \ + bbP_transform = new cTransform(params[5], &geo.C3_O3_P_angle, geo.R_O3_P, params_grad[5]); \ + this->groups.push_back(bbP); \ + this->transforms.push_back(bbP_transform); \ + nP = addNode(parentC, groups.back(), transforms.back()); + +#define ADD_O5 \ + bbO5 = makeAtom("O5'", firstAtomIndex+3, residueName, residueIndex, atoms_global); \ + bbO5_transform = new cTransform(params[0], &geo.O3_P_O5_angle, geo.R_P_O5, params_grad[0]); \ + this->groups.push_back(bbO5); \ + this->transforms.push_back(bbO5_transform); \ + nO5 = addNode(nP, groups.back(), transforms.back()); + +#define ADD_O5_5Prime \ + bbO5 = makeAtom("O5'", firstAtomIndex, residueName, residueIndex, atoms_global); \ + if(parentC==NULL) \ + bbO5_transform = new cTransform(&zero_const, &zero_const, zero_const, NULL); \ + else \ + bbO5_transform = new cTransform(params[0], &geo.O3_P_O5_angle, geo.R_P_O5, params_grad[0]); \ + this->groups.push_back(bbO5); \ + this->transforms.push_back(bbO5_transform); \ + nO5 = addNode(parentC, groups.back(), transforms.back()); + +#define ADD_C5 \ + bbC5 = makeAtom("C5'", firstAtomIndex+4, residueName, residueIndex, atoms_global); \ + bbC5_transform = new cTransform(params[1], &geo.P_O5_C5_angle, geo.R_O5_C5, params_grad[1]); \ + this->groups.push_back(bbC5); \ + this->transforms.push_back(bbC5_transform); \ + nC5 = addNode(nO5, groups.back(), transforms.back()); + +#define ADD_C5_5Prime \ + bbC5 = makeAtom("C5'", firstAtomIndex+1, residueName, residueIndex, atoms_global); \ + bbC5_transform = new cTransform(params[1], &geo.P_O5_C5_angle, geo.R_O5_C5, params_grad[1]); \ + this->groups.push_back(bbC5); \ + this->transforms.push_back(bbC5_transform); \ + nC5 = addNode(nO5, groups.back(), transforms.back()); + +#define ADD_C4 \ + bbC4 = makeAtom("C4'", firstAtomIndex+5, residueName, residueIndex, atoms_global); \ + bbC4_transform = new cTransform(params[2], &geo.O5_C5_C4_angle, geo.R_C5_C4, params_grad[2]); \ + this->groups.push_back(bbC4); \ + this->transforms.push_back(bbC4_transform); \ + nC4 = addNode(nC5, groups.back(), transforms.back()); + +#define ADD_C4_5Prime \ + bbC4 = makeAtom("C4'", firstAtomIndex+2, residueName, residueIndex, atoms_global); \ + bbC4_transform = new cTransform(params[2], &geo.O5_C5_C4_angle, geo.R_C5_C4, params_grad[2]); \ + this->groups.push_back(bbC4); \ + this->transforms.push_back(bbC4_transform); \ + nC4 = addNode(nC5, groups.back(), transforms.back()); + +#define ADD_C4_DT \ + cTransform *C4_dummy_transform = new cTransform(params[6], &geo.C4_correction_angle, 0.0, params_grad[6]); \ + cRigidGroup *C4_dummy_group = new cRigidGroup(); \ + this->groups.push_back(C4_dummy_group); \ + this->transforms.push_back(C4_dummy_transform); \ + cNode *C4_dummy_node = addNode(nC4, groups.back(), transforms.back()); + +#define ADD_O4 \ + bbO4 = makeAtom("O4'", firstAtomIndex+6, residueName, residueIndex, atoms_global); \ + bbO4_transform = new cTransform(params[10], &geo.C3_C4_O4_angle, geo.R_C4_O4, params_grad[10]); \ + this->groups.push_back(bbO4); \ + this->transforms.push_back(bbO4_transform); \ + nO4 = addNode(C4_dummy_node, groups.back(), transforms.back()); + +#define ADD_O4_5Prime \ + bbO4 = makeAtom("O4'", firstAtomIndex+3, residueName, residueIndex, atoms_global); \ + bbO4_transform = new cTransform(params[10], &geo.C3_C4_O4_angle, geo.R_C4_O4, params_grad[10]); \ + this->groups.push_back(bbO4); \ + this->transforms.push_back(bbO4_transform); \ + nO4 = addNode(C4_dummy_node, groups.back(), transforms.back()); + +#define ADD_C3 \ + bbC3 = makeAtom("C3'", firstAtomIndex+7, residueName, residueIndex, atoms_global); \ + bbC3_transform = new cTransform(params[3], &geo.C5_C4_C3_angle, geo.R_C4_C3, params_grad[3]); \ + this->groups.push_back(bbC3); \ + this->transforms.push_back(bbC3_transform); \ + nC3 = addNode(nC4, groups.back(), transforms.back()); + +#define ADD_C3_5Prime \ + bbC3 = makeAtom("C3'", firstAtomIndex+4, residueName, residueIndex, atoms_global); \ + bbC3_transform = new cTransform(params[3], &geo.C5_C4_C3_angle, geo.R_C4_C3, params_grad[3]); \ + this->groups.push_back(bbC3); \ + this->transforms.push_back(bbC3_transform); \ + nC3 = addNode(nC4, groups.back(), transforms.back()); + +#define ADD_C3_DT \ + cTransform *C3_dummy_transform = new cTransform(params[9], &geo.C3_correction_angle, 0.0, params_grad[9]); \ + cRigidGroup *C3_dummy_group = new cRigidGroup(); \ + this->groups.push_back(C3_dummy_group); \ + this->transforms.push_back(C3_dummy_transform); \ + cNode *C3_dummy_node = addNode(nC3, groups.back(), transforms.back()); + +#define ADD_C2 \ + bbC2 = makeAtom("C2'", firstAtomIndex+9, residueName, residueIndex, atoms_global); \ + bbC2_transform = new cTransform(params[8], &geo.O3_C3_C2_angle, geo.R_C3_C2, params_grad[8]); \ + this->groups.push_back(bbC2); \ + this->transforms.push_back(bbC2_transform); \ + nC2 = addNode(C3_dummy_node, groups.back(), transforms.back()); + +#define ADD_C2_5Prime \ + bbC2 = makeAtom("C2'", firstAtomIndex+6, residueName, residueIndex, atoms_global); \ + bbC2_transform = new cTransform(params[8], &geo.O3_C3_C2_angle, geo.R_C3_C2, params_grad[8]); \ + this->groups.push_back(bbC2); \ + this->transforms.push_back(bbC2_transform); \ + nC2 = addNode(C3_dummy_node, groups.back(), transforms.back()); + +#define ADD_C2GROUP \ + bbC2 = makeC2Group(geo, firstAtomIndex+9, residueName, residueIndex, atoms_global); \ + bbC2_transform = new cTransform(params[8], &geo.C3_O3_P_angle, geo.R_O3_P, params_grad[8]); \ + this->groups.push_back(bbC2); \ + this->transforms.push_back(bbC2_transform); \ + nC2 = addNode(C3_dummy_node, groups.back(), transforms.back()); + +#define ADD_C2GROUP_5Prime \ + bbC2 = makeC2Group(geo, firstAtomIndex+6, residueName, residueIndex, atoms_global); \ + bbC2_transform = new cTransform(params[8], &geo.C3_O3_P_angle, geo.R_O3_P, params_grad[8]); \ + this->groups.push_back(bbC2); \ + this->transforms.push_back(bbC2_transform); \ + nC2 = addNode(C3_dummy_node, groups.back(), transforms.back()); + +#define ADD_C1 \ + bbC1 = makeAtom("C1'", firstAtomIndex+10, residueName, residueIndex, atoms_global); \ + bbC1_transform = new cTransform(params[7], &geo.C3_C2_C1_angle, geo.R_C2_C1, params_grad[7]); \ + this->groups.push_back(bbC1); \ + this->transforms.push_back(bbC1_transform); \ + nC1 = addNode(nC2, groups.back(), transforms.back()); + +#define ADD_C1_5Prime \ + bbC1 = makeAtom("C1'", firstAtomIndex+7, residueName, residueIndex, atoms_global); \ + bbC1_transform = new cTransform(params[7], &geo.C3_C2_C1_angle, geo.R_C2_C1, params_grad[7]); \ + this->groups.push_back(bbC1); \ + this->transforms.push_back(bbC1_transform); \ + nC1 = addNode(nC2, groups.back(), transforms.back()); + +#define ADD_O3 \ + bbO3 = makeAtom("O3'", firstAtomIndex+8, residueName, residueIndex, atoms_global); \ + bbO3_transform = new cTransform(params[4], &geo.C4_C3_O3_angle, geo.R_C3_O3, params_grad[4]); \ + this->groups.push_back(bbO3); \ + this->transforms.push_back(bbO3_transform); \ + nO3 = addNode(nC3, groups.back(), transforms.back()); + +#define ADD_O3_5Prime \ + bbO3 = makeAtom("O3'", firstAtomIndex+5, residueName, residueIndex, atoms_global); \ + bbO3_transform = new cTransform(params[4], &geo.C4_C3_O3_angle, geo.R_C3_O3, params_grad[4]); \ + this->groups.push_back(bbO3); \ + this->transforms.push_back(bbO3_transform); \ + nO3 = addNode(nC3, groups.back(), transforms.back()); + +#define ADD_C1_R \ + bbC1 = makeAtom("C1'", firstAtomIndex+11, residueName, residueIndex, atoms_global); \ + bbC1_transform = new cTransform(params[7], &geo.C3_C2_C1_angle, geo.R_C2_C1, params_grad[7]); \ + this->groups.push_back(bbC1); \ + this->transforms.push_back(bbC1_transform); \ + nC1 = addNode(nC2, groups.back(), transforms.back()); + +#define ADD_C1_R_5Prime \ + bbC1 = makeAtom("C1'", firstAtomIndex+8, residueName, residueIndex, atoms_global); \ + bbC1_transform = new cTransform(params[7], &geo.C3_C2_C1_angle, geo.R_C2_C1, params_grad[7]); \ + this->groups.push_back(bbC1); \ + this->transforms.push_back(bbC1_transform); \ + nC1 = addNode(nC2, groups.back(), transforms.back()); + template cNode *cConformation::addGly(cNode *parentC, std::vector params, std::vector params_grad, bool terminal){ cNode *nC, *nCA, *nN; @@ -545,5 +756,499 @@ template cNode *cConformation::addTrp(cNode *parentC, std: return nC; } +template cNode *cConformation::addDG(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + cNode *nP, *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nN9, *nO3; + cTransform *bbP_transform, *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbN9_transform, *bbO3_transform; + cRigidGroup *bbP, *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbN9, *bbO3; + + uint residueIndex, firstAtomIndex; + char residueName = 'G'; + + PARENT_CHECK_DNA + + ADD_PHOSPHATE + ADD_O5 + ADD_C5 + ADD_C4 + ADD_C4_DT + ADD_O4 + ADD_C3 + ADD_C3_DT + ADD_C2 + ADD_C1 + + bbN9 = makeGuaGroup(geo, firstAtomIndex+11, residueName, residueIndex, atoms_global); + bbN9_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Gua_N1, params_grad[11]); //c_c_n_angle and c_n + this->groups.push_back(bbN9); + this->transforms.push_back(bbN9_transform); + nN9 = addNode(nC1, groups.back(), transforms.back()); + + ADD_O3 + + return nO3; // +} + +template cNode *cConformation::addDA(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + cNode *nP, *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nN9, *nO3; + cTransform *bbP_transform, *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbN9_transform, *bbO3_transform; + cRigidGroup *bbP, *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbN9, *bbO3; + + uint residueIndex, firstAtomIndex; + char residueName = 'A'; + + PARENT_CHECK_DNA + + ADD_PHOSPHATE + ADD_O5 + ADD_C5 + ADD_C4 + ADD_C4_DT + ADD_O4 + ADD_C3 + ADD_C3_DT + ADD_C2 + ADD_C1 + + bbN9 = makeAdeGroup(geo, firstAtomIndex+11, residueName, residueIndex, atoms_global); + bbN9_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Ade_N1, params_grad[11]); //c_c_n_angle and c_n + this->groups.push_back(bbN9); + this->transforms.push_back(bbN9_transform); + nN9 = addNode(nC1, groups.back(), transforms.back()); + + ADD_O3 + + return nO3; // +} + +template cNode *cConformation::addDT(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + cNode *nP, *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nO3,*nN1; + cTransform *bbP_transform, *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbO3_transform, *bbN1_transform; + cRigidGroup *bbP, *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbO3, *bbN1; + + uint residueIndex, firstAtomIndex; + char residueName = 'T'; + + PARENT_CHECK_DNA + + ADD_PHOSPHATE + ADD_O5 + ADD_C5 + ADD_C4 + ADD_C4_DT + ADD_O4 + ADD_C3 + ADD_C3_DT + ADD_C2 + ADD_C1 + ADD_O3 + + bbN1 = makeThyGroup(geo, firstAtomIndex+11, residueName, residueIndex, atoms_global); + bbN1_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Thy_N1, params_grad[11]); + this->transforms.push_back(bbN1_transform); + nN1 = addNode(nC1, groups.back(), transforms.back()); + + return nO3; // +} + +template cNode *cConformation::addDC(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + + cNode *nP, *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nO3, *nN1; + cTransform *bbP_transform, *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbO3_transform, *bbN1_transform; + cRigidGroup *bbP, *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbO3, *bbN1; + + uint residueIndex, firstAtomIndex; + char residueName = 'C'; + + PARENT_CHECK_DNA + + ADD_PHOSPHATE + ADD_O5 + ADD_C5 + ADD_C4 + ADD_C4_DT + ADD_O4 + ADD_C3 + ADD_C3_DT + ADD_C2 + ADD_C1 + ADD_O3 + + bbN1 = makeCytGroup(geo, firstAtomIndex+11, residueName, residueIndex, atoms_global); + bbN1_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Cyt_N1, params_grad[11]); + this->groups.push_back(bbN1); + this->transforms.push_back(bbN1_transform); + nN1 = addNode(nC1, groups.back(), transforms.back()); + + return nO3; // +} +template cNode *cConformation::addDG_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + cNode *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nN9, *nO3; + cTransform *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbN9_transform, *bbO3_transform; + cRigidGroup *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbN9, *bbO3; + + uint residueIndex, firstAtomIndex; + char residueName = 'G'; + + PARENT_CHECK_DNA + + ADD_O5_5Prime + ADD_C5_5Prime + ADD_C4_5Prime + ADD_C4_DT + ADD_O4_5Prime + ADD_C3_5Prime + ADD_C3_DT + ADD_C2_5Prime + ADD_C1_5Prime + ADD_O3_5Prime + + bbN9 = makeGuaGroup(geo, firstAtomIndex+8, residueName, residueIndex, atoms_global); + bbN9_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Gua_N1, params_grad[11]); //c_c_n_angle and c_n + this->groups.push_back(bbN9); + this->transforms.push_back(bbN9_transform); + nN9 = addNode(nC1, groups.back(), transforms.back()); + + return nO3; +} + +template cNode *cConformation::addDA_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + cNode *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nN9, *nO3; + cTransform *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbN9_transform, *bbO3_transform; + cRigidGroup *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbN9, *bbO3; + + uint residueIndex, firstAtomIndex; + char residueName = 'A'; + + PARENT_CHECK_DNA + + ADD_O5_5Prime + ADD_C5_5Prime + ADD_C4_5Prime + ADD_C4_DT + ADD_O4_5Prime + ADD_C3_5Prime + ADD_C3_DT + ADD_C2_5Prime + ADD_C1_5Prime + ADD_O3_5Prime + bbN9 = makeAdeGroup(geo, firstAtomIndex+8, residueName, residueIndex, atoms_global); + bbN9_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Ade_N1, params_grad[11]); //c_c_n_angle and c_n + this->groups.push_back(bbN9); + this->transforms.push_back(bbN9_transform); + nN9 = addNode(nC1, groups.back(), transforms.back()); + + return nO3; +} + +template cNode *cConformation::addDT_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + cNode *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nN1, *nO3; + cTransform *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbN1_transform, *bbO3_transform; + cRigidGroup *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbN1, *bbO3; + + uint residueIndex, firstAtomIndex; + char residueName = 'T'; + + PARENT_CHECK_DNA + + ADD_O5_5Prime + ADD_C5_5Prime + ADD_C4_5Prime + ADD_C4_DT + ADD_O4_5Prime + ADD_C3_5Prime + ADD_C3_DT + ADD_C2_5Prime + ADD_C1_5Prime + ADD_O3_5Prime + + bbN1 = makeThyGroup(geo, firstAtomIndex+8, residueName, residueIndex, atoms_global); + bbN1_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Thy_N1, params_grad[11]); //c_c_n_angle and c_n + this->groups.push_back(bbN1); + this->transforms.push_back(bbN1_transform); + nN1 = addNode(nC1, groups.back(), transforms.back()); + + return nO3; +} + +template cNode *cConformation::addDC_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + cNode *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nO3, *nN1; + cTransform *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbO3_transform, *bbN1_transform; + cRigidGroup *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbO3, *bbN1; + + uint residueIndex, firstAtomIndex; + char residueName = 'C'; + + PARENT_CHECK_DNA + + ADD_O5_5Prime + ADD_C5_5Prime + ADD_C4_5Prime + ADD_C4_DT + ADD_O4_5Prime + ADD_C3_5Prime + ADD_C3_DT + ADD_C2_5Prime + ADD_C1_5Prime + ADD_O3_5Prime + + bbN1 = makeCytGroup(geo, firstAtomIndex+8, residueName, residueIndex, atoms_global); + bbN1_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Cyt_N1, params_grad[11]); + this->groups.push_back(bbN1); + this->transforms.push_back(bbN1_transform); + nN1 = addNode(nC1, groups.back(), transforms.back()); + + return nO3; +} +template cNode *cConformation::addG(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + cNode *nP, *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nN9, *nO3; + cTransform *bbP_transform, *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbN9_transform, *bbO3_transform; + cRigidGroup *bbP, *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbN9, *bbO3; + + uint residueIndex, firstAtomIndex; + char residueName = 'G'; + + PARENT_CHECK_RNA + + ADD_PHOSPHATE + ADD_O5 + ADD_C5 + ADD_C4 + ADD_C4_DT + ADD_O4 + ADD_C3 + ADD_C3_DT + ADD_C2GROUP + ADD_C1_R + + bbN9 = makeGuaGroup(geo, firstAtomIndex+12, residueName, residueIndex, atoms_global); + bbN9_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Gua_N1, params_grad[11]); //c_c_n_angle and c_n for RNA + this->groups.push_back(bbN9); + this->transforms.push_back(bbN9_transform); + nN9 = addNode(nC1, groups.back(), transforms.back()); + + ADD_O3 + + return nO3; +} + +template cNode *cConformation::addA(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + cNode *nP, *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nN9, *nO3; + cTransform *bbP_transform, *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbN9_transform, *bbO3_transform; + cRigidGroup *bbP, *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbN9, *bbO3; + + uint residueIndex, firstAtomIndex; + char residueName = 'A'; + + PARENT_CHECK_RNA //? + + ADD_PHOSPHATE + ADD_O5 + ADD_C5 + ADD_C4 + ADD_C4_DT + ADD_O4 + ADD_C3 + ADD_C3_DT + ADD_C2GROUP + ADD_C1_R + + bbN9 = makeAdeGroup(geo, firstAtomIndex+12, residueName, residueIndex, atoms_global); + bbN9_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Ade_N1, params_grad[11]); //c_c_n_angle and c_n + this->groups.push_back(bbN9); + this->transforms.push_back(bbN9_transform); + nN9 = addNode(nC1, groups.back(), transforms.back()); + + ADD_O3 + + return nO3; // +} + +template cNode *cConformation::addU(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + cNode *nP, *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nO3,*nN1; + cTransform *bbP_transform, *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbO3_transform, *bbN1_transform; + cRigidGroup *bbP, *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbO3, *bbN1; + + uint residueIndex, firstAtomIndex; + char residueName = 'U'; + + PARENT_CHECK_RNA + + ADD_PHOSPHATE + ADD_O5 + ADD_C5 + ADD_C4 + ADD_C4_DT + ADD_O4 + ADD_C3 + ADD_C3_DT + ADD_C2GROUP + ADD_C1_R + + bbN1 = makeUraGroup(geo, firstAtomIndex+12, residueName, residueIndex, atoms_global); + bbN1_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Ura_N1, params_grad[11]); //c_c_n_angle and c_n + this->groups.push_back(bbN1); + this->transforms.push_back(bbN1_transform); + nN1 = addNode(nC1, groups.back(), transforms.back()); + + ADD_O3 + + return nO3; +} + +template cNode *cConformation::addC(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + + cNode *nP, *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nO3, *nN1; + cTransform *bbP_transform, *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbO3_transform, *bbN1_transform; + cRigidGroup *bbP, *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbO3, *bbN1; + + uint residueIndex, firstAtomIndex; + char residueName = 'C'; + + PARENT_CHECK_RNA + + ADD_PHOSPHATE + ADD_O5 + ADD_C5 + ADD_C4 + ADD_C4_DT + ADD_O4 + ADD_C3 + ADD_C3_DT + ADD_C2GROUP + ADD_C1_R + + bbN1 = makeCytGroup(geo, firstAtomIndex+12, residueName, residueIndex, atoms_global); + bbN1_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Cyt_N1, params_grad[11]); //c_c_n_angle and c_n + this->groups.push_back(bbN1); + this->transforms.push_back(bbN1_transform); + nN1 = addNode(nC1, groups.back(), transforms.back()); + + ADD_O3 + + return nO3; +} +template cNode *cConformation::addG_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + cNode *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nN9, *nO3; + cTransform *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbN9_transform, *bbO3_transform; + cRigidGroup *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbN9, *bbO3; + + uint residueIndex, firstAtomIndex; + char residueName = 'G'; + + PARENT_CHECK_RNA + + ADD_O5_5Prime + ADD_C5_5Prime + ADD_C4_5Prime + ADD_C4_DT + ADD_O4_5Prime + ADD_C3_5Prime + ADD_C3_DT + ADD_C2GROUP_5Prime + ADD_C1_R_5Prime + + bbN9 = makeGuaGroup(geo, firstAtomIndex+9, residueName, residueIndex, atoms_global); + bbN9_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Gua_N1, params_grad[11]); //c_c_n_angle and c_n + this->groups.push_back(bbN9); + this->transforms.push_back(bbN9_transform); + nN9 = addNode(nC1, groups.back(), transforms.back()); + + ADD_O3_5Prime + + return nO3; +} + +template cNode *cConformation::addA_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + cNode *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nN9, *nO3; + cTransform *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbN9_transform, *bbO3_transform; + cRigidGroup *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbN9, *bbO3; + + uint residueIndex, firstAtomIndex; + char residueName = 'A'; + + PARENT_CHECK_RNA + + ADD_O5_5Prime + ADD_C5_5Prime + ADD_C4_5Prime + ADD_C4_DT + ADD_O4_5Prime + ADD_C3_5Prime + ADD_C3_DT + ADD_C2GROUP_5Prime + ADD_C1_R_5Prime + + bbN9 = makeAdeGroup(geo, firstAtomIndex+9, residueName, residueIndex, atoms_global); + bbN9_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Ade_N1, params_grad[11]); //c_c_n_angle and c_n + this->groups.push_back(bbN9); + this->transforms.push_back(bbN9_transform); + nN9 = addNode(nC1, groups.back(), transforms.back()); + + ADD_O3_5Prime + + return nO3; +} + +template cNode *cConformation::addU_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + cNode *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nN1, *nO3; + cTransform *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbN1_transform, *bbO3_transform; + cRigidGroup *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbN1, *bbO3; + + uint residueIndex, firstAtomIndex; + char residueName = 'U'; + + PARENT_CHECK_RNA + + ADD_O5_5Prime + ADD_C5_5Prime + ADD_C4_5Prime + ADD_C4_DT + ADD_O4_5Prime + ADD_C3_5Prime + ADD_C3_DT + ADD_C2GROUP_5Prime + ADD_C1_R_5Prime + + bbN1 = makeUraGroup(geo, firstAtomIndex+9, residueName, residueIndex, atoms_global); + bbN1_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Ura_N1, params_grad[11]); //c_c_n_angle and c_n + this->groups.push_back(bbN1); + this->transforms.push_back(bbN1_transform); + nN1 = addNode(nC1, groups.back(), transforms.back()); + + ADD_O3_5Prime + + return nO3; +} + +template cNode *cConformation::addC_5Prime(cNode *parentC, std::vector params, std::vector params_grad, char last_res){ + cNode *nO5, *nC5, *nC4, *nO4, *nC3, *nC2, *nC1, *nO3, *nN1; + cTransform *bbO5_transform, *bbC5_transform, *bbC4_transform, *bbO4_transform, *bbC3_transform, *bbC2_transform, *bbC1_transform, *bbO3_transform, *bbN1_transform; + cRigidGroup *bbO5, *bbC5, *bbC4, *bbO4, *bbC3, *bbC2, *bbC1, *bbO3, *bbN1; + + uint residueIndex, firstAtomIndex; + char residueName = 'C'; + + PARENT_CHECK_RNA + + ADD_O5_5Prime + ADD_C5_5Prime + ADD_C4_5Prime + ADD_C4_DT + ADD_O4_5Prime + ADD_C3_5Prime + ADD_C3_DT + ADD_C2GROUP_5Prime + ADD_C1_R_5Prime + + bbN1 = makeCytGroup(geo, firstAtomIndex+9, residueName, residueIndex, atoms_global); + bbN1_transform = new cTransform(params[11], &geo.DNA_N1_angle, geo.R_Cyt_N1, params_grad[11]); //c_c_n_angle and c_n + this->groups.push_back(bbN1); + this->transforms.push_back(bbN1_transform); + nN1 = addNode(nC1, groups.back(), transforms.back()); + + ADD_O3_5Prime + + return nO3; +} template class cConformation; template class cConformation; \ No newline at end of file diff --git a/Layers/FullAtomModel/cGeometry.cpp b/Layers/FullAtomModel/cGeometry.cpp index f04090d..63de86a 100644 --- a/Layers/FullAtomModel/cGeometry.cpp +++ b/Layers/FullAtomModel/cGeometry.cpp @@ -5,84 +5,201 @@ #define KAPPA2 (3.14159 - 2.061) #define KAPPA3 (3.14159 -2.1186) #define OMEGACIS -3.1318 - template cGeometry::cGeometry(){ - //backbone angles - C_N_CA_angle = (M_PI - 1.9391); - N_CA_C_angle = (M_PI - 2.061); - CA_C_N_angle = (M_PI - 2.1186); - CA_C_O_angle = (M_PI - 2.1033); + //backbone angles [Original TPL angles: from Peptide Builder] +// C_N_CA_angle = (M_PI - 1.9391); +// N_CA_C_angle = (M_PI - 2.061); +// CA_C_N_angle = (M_PI - 2.1186); +// CA_C_O_angle = (M_PI - 2.1033); +// omega_const = -3.1318; + //backbone angles [Engh & Huber angles] + C_N_CA_angle = (M_PI - 2.1241); + N_CA_C_angle = (M_PI - 2.0281); + CA_C_N_angle = (M_PI - 1.9408); + CA_C_O_angle = (M_PI - 2.0420); omega_const = -3.1318; + //backbone distatnces +// R_CA_C = 1.525; +// R_C_N = 1.330; +// R_N_CA = 1.460; + +//backbone distatnces [Engh & Huber] R_CA_C = 1.525; - R_C_N = 1.330; - R_N_CA = 1.460; + R_C_N = 1.329; + R_N_CA = 1.458; + //C-beta - R_CA_CB = 1.52; +// R_CA_CB = 1.52; + // C_CA_CB_angle = -(M_PI - 1.9111); +// C_CA_CB_angle = -(M_PI -M_PI*109.5/180.0); +// N_C_CA_CB_diangle = 0.5*M_PI*122.6860/180.0; +// correction_angle = 0.0;//-C_CA_CB_angle/2.0; + + //C-beta + R_CA_CB = 1.53; // C_CA_CB_angle = -(M_PI - 1.9111); - C_CA_CB_angle = -(M_PI -M_PI*109.5/180.0); + C_CA_CB_angle = -(M_PI -M_PI*110.1/180.0); N_C_CA_CB_diangle = 0.5*M_PI*122.6860/180.0; correction_angle = 0.0;//-C_CA_CB_angle/2.0; +//Original TPL values //OG serine - R_CB_OG = 1.417; - CA_CB_OG_angle = M_PI*110.773/180.0; +// R_CB_OG = 1.417; +// CA_CB_OG_angle = M_PI*110.773/180.0; //SG cysteine - R_CB_SG = 1.808; - CA_CB_SG_angle = M_PI*113.8169/180.0; +// R_CB_SG = 1.808; +// CA_CB_SG_angle = M_PI*113.8169/180.0; //Valine CG1 and CG2 - R_CB_CG = 1.527; - CA_CB_CG1_angle = M_PI*110.7/180.0; - CG1_CB_CG2_angle = 1.88444687881; +// R_CB_CG = 1.527; +// CA_CB_CG1_angle = M_PI*110.7/180.0; +// CG1_CB_CG2_angle = 1.88444687881; //Isoleucine - R_CG1_CD1 = 1.52; - CB_CG1_CD1_angle = M_PI*113.97/180.0; +// R_CG1_CD1 = 1.52; +// CB_CG1_CD1_angle = M_PI*113.97/180.0; //Threonine - R_CB_OG1 = 1.43; - CA_CB_OG1_angle = M_PI*109.18/180.0; - OG1_CB_CG2_angle = 1.90291866134; +// R_CB_OG1 = 1.43; +// CA_CB_OG1_angle = M_PI*109.18/180.0; +// OG1_CB_CG2_angle = 1.90291866134; //Arginine +// R_C_C = R_CA_CB; +// C_C_C_angle = C_CA_CB_angle; +// R_CD_NE = 1.46; +// R_NE_CZ = 1.33; +// R_CZ_NH = 1.33; +// CA_CB_CG_angle = -(M_PI - M_PI*113.83/180.0); +// CB_CG_CD_angle = -(M_PI - M_PI*111.79/180.0); +// CG_CD_NE_angle = -(M_PI - M_PI*111.68/180.0); +// CD_NE_CZ_angle = -(M_PI - M_PI*124.79/180.0); +// NE_CZ_NH_angle = M_PI*124.79/180.0; + //Lysine +// R_CD_CE = 1.46; +// R_CE_NZ = 1.33; +// CD_CE_NZ_angle = -(M_PI - M_PI*124.79/180.0); + //Aspartic acid +// R_CG_OD = 1.25; +// CB_CG_OD_angle = M_PI*119.22/180.0; +// OD1_CG_OD2_angle = 2.13911043783; + //Asparagine +// R_CG_OD1 = 1.23; +// R_CG_ND2 = 1.33; +// CB_CG_OD1_angle = M_PI*120.85/180.0; +// CB_CG_ND2_angle = M_PI*116.48/180.0; + //Methionine +// R_CG_SD = 1.81; +// R_SD_CE = 1.79; +// CB_CG_SD_angle = -(M_PI - M_PI*112.69/180.0); +// CG_SD_CE_angle = -(M_PI - M_PI*100.61/180.0); +// //Histidine +// R_CG_CD2 = 1.35; +// R_CG_ND1 = 1.38; +// R_CG_NE2 = 2.1912; +// R_CG_CE1 = 2.1915; +// CB_CG_CD2_angle = -(M_PI - 2.27957453603); +// CB_CG_ND1_angle = (M_PI - 2.14413698608); +// CB_CG_NE2_angle = -(M_PI - 2.90352974362); +// CB_CG_CE1_angle = (M_PI - 2.75209584734); + + +// Engh & Huber + //OG serine + R_CB_OG = 1.417; + CA_CB_OG_angle = M_PI*111.1/180.0; + //SG cysteine + R_CB_SG = 1.803; + CA_CB_SG_angle = M_PI*114.4/180.0; + //Valine CG1 and CG2 + R_CB_CG = 1.521; + CA_CB_CG1_angle = M_PI*110.5/180.0; + CG1_CB_CG2_angle = 1.93382481121; + //Isoleucine + R_CG1_CD1 = 1.513; + CB_CG1_CD1_angle = M_PI*113.8/180.0; + //Threonine + R_CB_OG1 = 1.433; + CA_CB_OG1_angle = M_PI*109.6/180.0; + OG1_CB_CG2_angle = 1.90764487243; + //Arginine R_C_C = R_CA_CB; C_C_C_angle = C_CA_CB_angle; - R_CD_NE = 1.46; - R_NE_CZ = 1.33; - R_CZ_NH = 1.33; - CA_CB_CG_angle = -(M_PI - M_PI*113.83/180.0); - CB_CG_CD_angle = -(M_PI - M_PI*111.79/180.0); - CG_CD_NE_angle = -(M_PI - M_PI*111.68/180.0); - CD_NE_CZ_angle = -(M_PI - M_PI*124.79/180.0); - NE_CZ_NH_angle = M_PI*124.79/180.0; + R_CD_NE = 1.460; + R_NE_CZ = 1.329; + R_CZ_NH = 1.326; + CA_CB_CG_angle = -(M_PI - M_PI*111.2/180.0); + CB_CG_CD_angle = -(M_PI - M_PI*111.3/180.0); + CG_CD_NE_angle = -(M_PI - M_PI*112.0/180.0); + CD_NE_CZ_angle = -(M_PI - M_PI*124.2/180.0); + NE_CZ_NH_angle = M_PI*120.0/180.0; //Lysine - R_CD_CE = 1.46; - R_CE_NZ = 1.33; - CD_CE_NZ_angle = -(M_PI - M_PI*124.79/180.0); + R_CD_CE = 1.52; + R_CE_NZ = 1.326; + CD_CE_NZ_angle = -(M_PI - M_PI*111.9/180.0); //Aspartic acid - R_CG_OD = 1.25; - CB_CG_OD_angle = M_PI*119.22/180.0; - OD1_CG_OD2_angle = 2.13911043783; + R_CG_OD = 1.249; + CB_CG_OD_angle = M_PI*118.4/180.0; + OD1_CG_OD2_angle = 2.14500965070; //Asparagine - R_CG_OD1 = 1.23; - R_CG_ND2 = 1.33; + R_CG_OD1 = 1.231; + R_CG_ND2 = 1.328; CB_CG_OD1_angle = M_PI*120.85/180.0; - CB_CG_ND2_angle = M_PI*116.48/180.0; + CB_CG_ND2_angle = M_PI*116.4/180.0; //Methionine - R_CG_SD = 1.81; - R_SD_CE = 1.79; - CB_CG_SD_angle = -(M_PI - M_PI*112.69/180.0); - CG_SD_CE_angle = -(M_PI - M_PI*100.61/180.0); + R_CG_SD = 1.803; + R_SD_CE = 1.791; + CB_CG_SD_angle = -(M_PI - M_PI*112.7/180.0); + CG_SD_CE_angle = -(M_PI - M_PI*100.9/180.0); //Histidine - R_CG_CD2 = 1.35; - R_CG_ND1 = 1.38; - R_CG_NE2 = 2.1912; - R_CG_CE1 = 2.1915; - CB_CG_CD2_angle = -(M_PI - 2.27957453603); - CB_CG_ND1_angle = (M_PI - 2.14413698608); - CB_CG_NE2_angle = -(M_PI - 2.90352974362); - CB_CG_CE1_angle = (M_PI - 2.75209584734); - - - + R_CG_CD2 = 1.356; + R_CG_ND1 = 1.371; + R_CG_NE2 = 2.187; + R_CG_CE1 = 2.167;//avg from 2.14289 + 2.19018 + CB_CG_CD2_angle = -(M_PI - 2.25322006432); + CB_CG_ND1_angle = (M_PI - 2.12232037043); + CB_CG_NE2_angle = -(M_PI - 2.89899188756); + CB_CG_CE1_angle = (M_PI - 2.79078147394); + + + //Geometry for Nucleic Acids + //NA Backbone + R_P_O5 = 1.593; + R_O5_C5 = 1.440; + R_C5_C4 = 1.511; // poly_type + R_C4_C3 = 1.528; // poly_type + R_C3_O3 = 1.431; // poly_type + R_O3_P = 1.607; + O3_P_O5_angle = (M_PI - (M_PI*104.0/180)); + P_O5_C5_angle = (M_PI - (M_PI*120.9/180)); + O5_C5_C4_angle = (M_PI - (M_PI*110.2/180)); // poly_type + C5_C4_C3_angle = (M_PI - (M_PI*114.7/180)); // poly_type + C4_C3_O3_angle = (M_PI - (M_PI*110.3/180)); // poly_type + C3_O3_P_angle = (M_PI - (M_PI*119.7/180)); + + //NA sugar ring + R_C4_O4 = 1.446; // poly_type + R_C3_C2 = 1.525; // poly_type + R_C2_C1 = 1.528; // poly_type + + C3_C4_O4_angle = (M_PI - (M_PI*109.4/180)); //poly_type C5_C4_O4 109.4 || C3_C4_O4 105.6 + O3_C3_C2_angle = (M_PI - (M_PI*103.2/180)); //poly_type C4_C3_C2 103.2 || O3_C3_C2 110.6 + C3_C2_C1_angle = (M_PI - (M_PI*102.7/180)); //poly_type + + //C4_dummy_transform + C4_correction_angle = 0.0; + O3_C3_C4_O4_diangle = -0.5*((M_PI*252.20807/180)); //approximation from Gaussian + + //C3_dummy_transform + C3_correction_angle = 0.0; + P_O3_C3_C2_diangle = -0.5*(M_PI*93.58720/180); //approximation from Gaussian + +// DNA chi + R_Cyt_N1 = 1.470; + R_Gua_N1 = 1.459; + R_Thy_N1 = 1.473; + R_Ade_N1 = 1.462; + R_Ura_N1 = 1.469; + DNA_N1_angle = (M_PI - (M_PI * 114.2/180)); + } template cGeometry::~cGeometry(){ @@ -90,34 +207,66 @@ template cGeometry::~cGeometry(){ } template void cGeometry::gly(){ - //backbone angles - C_N_CA_angle = (3.14159 - 1.9391); - N_CA_C_angle = (3.14159 - 2.061); - CA_C_N_angle = (3.14159 - 2.1186); - CA_C_O_angle = (3.14159 - 2.1033); + //backbone angles [Original TPL values] +// C_N_CA_angle = (3.14159 - 1.9391); +// N_CA_C_angle = (3.14159 - 2.061); +// CA_C_N_angle = (3.14159 - 2.116); +// CA_C_O_angle = (3.14159 - 2.1033); +// omega_const = -3.1318; + +//backbone distatnces +// R_CA_C = 1.525; +// R_C_N = 1.330; +// R_N_CA = 1.460; + + + //backbone angles [Engh & Huber angles] + C_N_CA_angle = (M_PI - 2.1049); + N_CA_C_angle = (M_PI - 1.9635); + CA_C_N_angle = (M_PI - 2.0316); + CA_C_O_angle = (M_PI - 2.0665); omega_const = -3.1318; //backbone distatnces - R_CA_C = 1.525; - R_C_N = 1.330; - R_N_CA = 1.460; + R_CA_C = 1.516; + R_C_N = 1.329; + R_N_CA = 1.451; } template void cGeometry::ala(){ - //backbone angles - C_N_CA_angle = (3.14159 - 1.9391); - N_CA_C_angle = (3.14159 - 2.061); - CA_C_N_angle = (3.14159 - 2.1186); - CA_C_O_angle = (3.14159 - 2.1033); + //backbone angles [Original TPL values] +// C_N_CA_angle = (3.14159 - 1.9391); +// N_CA_C_angle = (3.14159 - 2.061); +// CA_C_N_angle = (3.14159 - 2.1186); +// CA_C_O_angle = (3.14159 - 2.1033); +// omega_const = -3.1318; + + //backbone angles [Engh & Huber angles] + C_N_CA_angle = (M_PI - 2.1241); + N_CA_C_angle = (M_PI - 2.0281); + CA_C_N_angle = (M_PI - 1.9408); + CA_C_O_angle = (M_PI - 2.0420); omega_const = -3.1318; - //backbone distatnces + + //backbone distatnces [Orig TPL] +// R_CA_C = 1.525; +// R_C_N = 1.330; +// R_N_CA = 1.460; + + //backbone distatnces [E&H] R_CA_C = 1.525; - R_C_N = 1.330; - R_N_CA = 1.460; + R_C_N = 1.329; + R_N_CA = 1.458; - //C-beta - R_CA_CB = 1.52; - C_CA_CB_angle = -(3.14159 - 1.9111);//(3.14159 - 1.9111); - N_C_CA_CB_diangle = M_PI*111.0/180.0;; + + //C-beta [Original TPL values] +// R_CA_CB = 1.52; +// C_CA_CB_angle = -(3.14159 - 1.9111);//(3.14159 - 1.9111); +// N_C_CA_CB_diangle = M_PI*111.0/180.0;; + + //C-beta [Engh and Huber Angles] + R_CA_CB = 1.521; + C_CA_CB_angle = -(3.14159 - 1.9286);//(3.14159 - 1.9111); + N_C_CA_CB_diangle = M_PI*111.0/180.0;; // need dihedral angle } template class cGeometry; diff --git a/Layers/FullAtomModel/cGeometry.h b/Layers/FullAtomModel/cGeometry.h index be7056c..cc7e014 100644 --- a/Layers/FullAtomModel/cGeometry.h +++ b/Layers/FullAtomModel/cGeometry.h @@ -52,6 +52,21 @@ class cGeometry{ T R_CG_CD2, R_CG_ND1, R_CG_NE2, R_CG_CE1; T CB_CG_CD2_angle, CB_CG_ND1_angle, CB_CG_NE2_angle, CB_CG_CE1_angle; + //NA Geometry + //NA backbone + T R_P_O5, R_O5_C5, R_C5_C4, R_C4_C3, R_C3_O3, R_O3_P; + T O3_P_O5_angle, P_O5_C5_angle, O5_C5_C4_angle, C5_C4_C3_angle, C4_C3_O3_angle, C3_O3_P_angle; + + //NA sugar ring + T R_C4_O4, R_C3_C2, R_C2_C1; + T C3_C4_O4_angle, O3_C3_C2_angle, C3_C2_C1_angle; + + // Sugar dummy transforms + T C4_correction_angle, O3_C3_C4_O4_diangle; + T C3_correction_angle, P_O3_C3_C2_diangle; + + T R_Cyt_N1, R_Gua_N1, R_Thy_N1, R_Ade_N1, R_Ura_N1, DNA_N1_angle; + cGeometry(); ~cGeometry(); diff --git a/Layers/FullAtomModel/cPDBLoader.cpp b/Layers/FullAtomModel/cPDBLoader.cpp index d813236..613625f 100644 --- a/Layers/FullAtomModel/cPDBLoader.cpp +++ b/Layers/FullAtomModel/cPDBLoader.cpp @@ -22,34 +22,127 @@ using namespace StringUtil; cPDBLoader::cPDBLoader(){ } -cPDBLoader::cPDBLoader(std::string filename) { +cPDBLoader::cPDBLoader(std::string filename, std::string chain_id, int polymer_type) { + int model_count; std::ifstream pfile(filename); std::string line, header, xStr, yStr, zStr, atom_name, res_name, chain_name; int res_num; // reading raw file + while ( getline (pfile,line) ){ header = line.substr(0,4); - - if( header.compare("ATOM")==0){ - atom_name = trim(line.substr(12,4)); - // std::cout<(std::stof(xStr),std::stof(yStr),std::stof(zStr))); - chain_names.push_back(chain_name); - res_names.push_back(res_name); - res_nums.push_back(res_num); - atom_names.push_back(atom_name); - // std::cout<(std::stof(xStr),std::stof(yStr),std::stof(zStr))); + chain_names.push_back(chain_name); + res_names.push_back(res_name); + res_nums.push_back(res_num); + atom_names.push_back(atom_name); + } + if(chain_name == chain_id && chain_id.length() > 0){ + xStr = line.substr(30,8); + yStr = line.substr(38,8); + zStr = line.substr(46,8); + res_name = trim(line.substr(17,3)); + res_num = std::stoi(line.substr(22,4)); + r.push_back(cVector3(std::stof(xStr),std::stof(yStr),std::stof(zStr))); + chain_names.push_back(chain_name); + res_names.push_back(res_name); + res_nums.push_back(res_num); + atom_names.push_back(atom_name); + } + } + } + + } + + if(polymer_type == 1 && model_count == 0){ + if( header.compare("ATOM")==0){ + atom_name = trim(line.substr(12,4)); + if(isHeavyAtom(atom_name)){ + res_name = trim(line.substr(17,3)); + if(isNucleotide(res_name, polymer_type)){ + chain_name = line.substr(21, 1); + if(chain_id.empty()){ + xStr = line.substr(30,8); + yStr = line.substr(38,8); + zStr = line.substr(46,8); + res_num = std::stoi(line.substr(22,4)); + r.push_back(cVector3(std::stof(xStr),std::stof(yStr),std::stof(zStr))); + chain_names.push_back(chain_name); + res_names.push_back(res_name); + res_nums.push_back(res_num); + atom_names.push_back(atom_name); + } + if(chain_name == chain_id && chain_id.length() > 0){ + xStr = line.substr(30,8); + yStr = line.substr(38,8); + zStr = line.substr(46,8); + res_num = std::stoi(line.substr(22,4)); + r.push_back(cVector3(std::stof(xStr),std::stof(yStr),std::stof(zStr))); + chain_names.push_back(chain_name); + res_names.push_back(res_name); + res_nums.push_back(res_num); + atom_names.push_back(atom_name); + } + } + } + } + } + if(polymer_type == 2 && model_count == 0){ + if( header.compare("ATOM")==0){ + atom_name = trim(line.substr(12,4)); + if(isHeavyAtom(atom_name)){ + res_name = trim(line.substr(17,3)); + if(isNucleotide(res_name, polymer_type)){ + chain_name = line.substr(21, 1); + if(chain_id.empty()){ + xStr = line.substr(30,8); + yStr = line.substr(38,8); + zStr = line.substr(46,8); + res_num = std::stoi(line.substr(22,4)); + r.push_back(cVector3(std::stof(xStr),std::stof(yStr),std::stof(zStr))); + chain_names.push_back(chain_name); + res_names.push_back(res_name); + res_nums.push_back(res_num); + atom_names.push_back(atom_name); + } + if (chain_name == chain_id && chain_id.length() > 0){ + xStr = line.substr(30,8); + yStr = line.substr(38,8); + zStr = line.substr(46,8); + res_num = std::stoi(line.substr(22,4)); + r.push_back(cVector3(std::stof(xStr),std::stof(yStr),std::stof(zStr))); + chain_names.push_back(chain_name); + res_names.push_back(res_name); + res_nums.push_back(res_num); + atom_names.push_back(atom_name); + } + } + } + } + } } } @@ -57,58 +150,49 @@ cPDBLoader::~cPDBLoader() { } -void cPDBLoader::reorder(){ - int global_ind=0; - int local_ind; - std::string lastO("O"); - - // std::cout<<"==="< cPDBLoader::getCenterMass(){ cVector3 c(0., 0., 0.); diff --git a/Layers/FullAtomModel/cPDBLoader.h b/Layers/FullAtomModel/cPDBLoader.h index c849f5a..facb937 100644 --- a/Layers/FullAtomModel/cPDBLoader.h +++ b/Layers/FullAtomModel/cPDBLoader.h @@ -24,7 +24,7 @@ class cPDBLoader { cVector3 b0, b1; public: cPDBLoader(); - cPDBLoader(std::string filename); + cPDBLoader(std::string filename, std::string chain_id=std::string(), int polymer_type=0); virtual ~cPDBLoader(); //order according to cConformation diff --git a/Layers/FullAtomModel/cRigidGroup.cpp b/Layers/FullAtomModel/cRigidGroup.cpp index a99ee30..aee06c7 100644 --- a/Layers/FullAtomModel/cRigidGroup.cpp +++ b/Layers/FullAtomModel/cRigidGroup.cpp @@ -209,6 +209,93 @@ template cRigidGroup *makeTrpGroup(cGeometry &geo, uint atomI return g; } +template cRigidGroup *makePhosGroup(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr){ + cRigidGroup *g = new cRigidGroup(); + g->addAtom(cVector3(0,0,0), "P", atomIndex, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(0.56, -1.24, 0.64), "OP1", atomIndex+1, residueName, residueIndex, atoms_global_ptr); // was exp(-0.014, 1.375, -0.45) now from guassian -> OP1 + g->addAtom(cVector3(0.32, 1.24, 0.79), "OP2", atomIndex+2, residueName, residueIndex, atoms_global_ptr); // was OP1(-0.17653, -1.23469, -0.805979) -> OP2 [original OP2 was (0.0135, -0.09, -1.5)] now from gauss + return g; +} +// g->addAtom(cVector3(-0.184111, -1.287736, 0.84606), "OP1", atomIndex+1, residueName, residueIndex, atoms_global_ptr); // was exp(-0.014, 1.375, -0.45) now from guassian -> OP1 +// g->addAtom(cVector3(1.40141, 0.67988, 0.365802), "OP2", atomIndex+2, residueName, residueIndex, atoms_global_ptr); // was OP1(-0.17653, -1.23469, -0.805979) -> OP2 [original OP2 was (0.0135, -0.09, -1.5)] now from gauss +template cRigidGroup *makeC2Group(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr){ + cRigidGroup *g = new cRigidGroup(); + g->addAtom(cVector3(0,0,0), "C2'", atomIndex, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(0.56, -1.24, 0.64), "O2'", atomIndex+1, residueName, residueIndex, atoms_global_ptr);//(0, 1.413, 0) or (0, 0, 1.413)? + return g; +} + +template cRigidGroup *makeCytGroup(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr){ + cRigidGroup *g = new cRigidGroup(); + g->addAtom(cVector3(0,0,0), "N1", atomIndex, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(0.688, -1.216, 0), "C2", atomIndex+1, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(0.0405, -2.2735, 0), "O2", atomIndex+2, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(2.041, -1.209, 0), "N3", atomIndex+3, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(2.7035, -0.056, 0), "C4", atomIndex+4, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(4.0375, -0.1026, 0), "N4", atomIndex+5, residueName, residueIndex, atoms_global_ptr); //changed for test, was (3.8475, -0.056, 0) and previously was (3.8475, 0.632, 0) + g->addAtom(cVector3(2.027, 1.2, 0), "C5", atomIndex+6, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(0.688, 1.181, 0), "C6", atomIndex+7, residueName, residueIndex, atoms_global_ptr); + return g; +} + +template cRigidGroup *makeThyGroup(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr){ + cRigidGroup *g = new cRigidGroup(); + g->addAtom(cVector3(0,0,0), "N1", atomIndex, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(0.675, -1.199, 0), "C2", atomIndex+1, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(0.0555, -2.25, 0), "O2", atomIndex+2, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(2.043, -1.084, 0), "N3", atomIndex+3, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(2.783, -0.073, 0), "C4", atomIndex+4, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(4.0063, -0.03403, 0), "O4", atomIndex+5, residueName, residueIndex, atoms_global_ptr); //changed for test, was (3.835, -0.073, 0) and previously was (3.835, -0.5595, 0) + g->addAtom(cVector3(2.009, 1.305, 0), "C5", atomIndex+6, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(2.7205, 2.6205, 0), "C7", atomIndex+7, residueName, residueIndex, atoms_global_ptr); //changed for test, was (2.0277205, 2.6205, 0) + g->addAtom(cVector3(0.674, 1.202, 0), "C6", atomIndex+8, residueName, residueIndex, atoms_global_ptr); + return g; +} + +template cRigidGroup *makeAdeGroup(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr){ + cRigidGroup *g = new cRigidGroup(); + g->addAtom(cVector3(0,0,0), "N9", atomIndex, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(0.828, 1.095, 0), "C8", atomIndex+1, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(2.105, 0.793, 0), "N7", atomIndex+2, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(2.095, -0.595, 0), "C5", atomIndex+3, residueName, residueIndex, atoms_global_ptr); // changed for test, was (1.95, -0.2865, 0) + g->addAtom(cVector3(3.2, -1.55, 0), "C6", atomIndex+4, residueName, residueIndex, atoms_global_ptr);// test, was (2.515, -1.5735, 0) + g->addAtom(cVector3(4.455, -1.25, 0), "N6", atomIndex+5, residueName, residueIndex, atoms_global_ptr);// test, was (3.68, -2.2245, 0) + g->addAtom(cVector3(2.79, -2.84, 0), "N1", atomIndex+6, residueName, residueIndex, atoms_global_ptr); // test, was (1.803, -3.3725, 0) + g->addAtom(cVector3(1.5, -3.2, 0), "C2", atomIndex+7, residueName, residueIndex, atoms_global_ptr);//test, was (0.798, -3.5, 0) + g->addAtom(cVector3(0.397, -2.427, 0), "N3", atomIndex+8, residueName, residueIndex, atoms_global_ptr);//test, was (0.227, -2.298, 0) + g->addAtom(cVector3(0.829, -1.096, 0), "C4", atomIndex+9, residueName, residueIndex, atoms_global_ptr); + return g; +} + +template cRigidGroup *makeGuaGroup(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr){ + cRigidGroup *g = new cRigidGroup(); + g->addAtom(cVector3(0,0,0), "N9", atomIndex, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(0.824, 1.101, 0), "C8", atomIndex+1, residueName, residueIndex, atoms_global_ptr); //(0.824, 1.101, 0) + g->addAtom(cVector3(2.092, 0.792, 0), "N7", atomIndex+2, residueName, residueIndex, atoms_global_ptr); //(2.092, 0.792, 0) + g->addAtom(cVector3(2.092, -0.596, 0), "C5", atomIndex+3, residueName, residueIndex, atoms_global_ptr); //(1.921, -0.266, 0) + g->addAtom(cVector3(3.18, -1.5, 0), "C6", atomIndex+4, residueName, residueIndex, atoms_global_ptr);// test, was (2.489, -1.566, 0) + g->addAtom(cVector3(4.356, -1.23, 0), "O6", atomIndex+5, residueName, residueIndex, atoms_global_ptr);// test, was (3.585, -2.139, 0) + g->addAtom(cVector3(2.8, -2.8, 0), "N1", atomIndex+6, residueName, residueIndex, atoms_global_ptr);// test, was (1.699, -2.711, 0) + g->addAtom(cVector3(1.5, -3.2, 0), "C2", atomIndex+7, residueName, residueIndex, atoms_global_ptr);// test, was (0.756333, -3.399666, 0) + g->addAtom(cVector3(1.28, -4.5, 0), "N2", atomIndex+8, residueName, residueIndex, atoms_global_ptr);// test, was (0.0651666, -4.5481666, 0) + g->addAtom(cVector3(0.385, -2.38, 0), "N3", atomIndex+9, residueName, residueIndex, atoms_global_ptr);//test, was (0.204, -2.3, 0) + g->addAtom(cVector3(0.823, -1.1, 0), "C4", atomIndex+10, residueName, residueIndex, atoms_global_ptr); + return g; +} + +template cRigidGroup *makeUraGroup(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr){ + cRigidGroup *g = new cRigidGroup(); + g->addAtom(cVector3(0,0,0), "N1", atomIndex, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(0.688, -1.216, 0), "C2", atomIndex+1, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(0.0405, -2.2735, 0), "O2", atomIndex+2, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(2.041, -1.209, 0), "N3", atomIndex+3, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(2.7035, -0.056, 0), "C4", atomIndex+4, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(4.0375, -0.1026, 0), "O4", atomIndex+5, residueName, residueIndex, atoms_global_ptr); //Will need to be adjusted + g->addAtom(cVector3(2.027, 1.2, 0), "C5", atomIndex+6, residueName, residueIndex, atoms_global_ptr); + g->addAtom(cVector3(0.688, 1.181, 0), "C6", atomIndex+7, residueName, residueIndex, atoms_global_ptr); + return g; +} + template class cRigidGroup; template cRigidGroup *makeAtom(std::string, uint, char, uint, float*); template cRigidGroup *makeCarbonyl(cGeometry&, uint, char, uint, float*, bool); @@ -227,6 +314,13 @@ template cRigidGroup *makeProGroup(cGeometry&, uint, char, uint, f template cRigidGroup *makePheGroup(cGeometry&, uint, char, uint, float*); template cRigidGroup *makeTyrGroup(cGeometry&, uint, char, uint, float*); template cRigidGroup *makeTrpGroup(cGeometry&, uint, char, uint, float*); +template cRigidGroup *makePhosGroup(cGeometry&, uint, char, uint, float*); +template cRigidGroup *makeC2Group(cGeometry&, uint, char, uint, float*); +template cRigidGroup *makeCytGroup(cGeometry&, uint, char, uint, float*); +template cRigidGroup *makeThyGroup(cGeometry&, uint, char, uint, float*); +template cRigidGroup *makeAdeGroup(cGeometry&, uint, char, uint, float*); +template cRigidGroup *makeGuaGroup(cGeometry&, uint, char, uint, float*); +template cRigidGroup *makeUraGroup(cGeometry&, uint, char, uint, float*); template class cRigidGroup; template cRigidGroup *makeAtom(std::string, uint, char, uint, double*); @@ -246,3 +340,10 @@ template cRigidGroup *makeProGroup(cGeometry&, uint, char, uint, template cRigidGroup *makePheGroup(cGeometry&, uint, char, uint, double*); template cRigidGroup *makeTyrGroup(cGeometry&, uint, char, uint, double*); template cRigidGroup *makeTrpGroup(cGeometry&, uint, char, uint, double*); +template cRigidGroup *makePhosGroup(cGeometry&, uint, char, uint, double*); +template cRigidGroup *makeC2Group(cGeometry&, uint, char, uint, double*); +template cRigidGroup *makeCytGroup(cGeometry&, uint, char, uint, double*); +template cRigidGroup *makeThyGroup(cGeometry&, uint, char, uint, double*); +template cRigidGroup *makeAdeGroup(cGeometry&, uint, char, uint, double*); +template cRigidGroup *makeGuaGroup(cGeometry&, uint, char, uint, double*); +template cRigidGroup *makeUraGroup(cGeometry&, uint, char, uint, double*); diff --git a/Layers/FullAtomModel/cRigidGroup.h b/Layers/FullAtomModel/cRigidGroup.h index 8c62426..948392e 100644 --- a/Layers/FullAtomModel/cRigidGroup.h +++ b/Layers/FullAtomModel/cRigidGroup.h @@ -37,6 +37,13 @@ template cRigidGroup *makeProGroup(cGeometry &geo, uint atomI template cRigidGroup *makePheGroup(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr); template cRigidGroup *makeTyrGroup(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr); template cRigidGroup *makeTrpGroup(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr); +template cRigidGroup *makePhosGroup(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr); +template cRigidGroup *makeC2Group(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr); +template cRigidGroup *makeCytGroup(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr); +template cRigidGroup *makeThyGroup(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr); +template cRigidGroup *makeAdeGroup(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr); +template cRigidGroup *makeGuaGroup(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr); +template cRigidGroup *makeUraGroup(cGeometry &geo, uint atomIndex, char residueName, uint residueIndex, T *atoms_global_ptr); template std::ostream& operator<<(std::ostream& os, const cRigidGroup& rg); diff --git a/Layers/FullAtomModel/main.cpp b/Layers/FullAtomModel/main.cpp index fe4df82..237cac3 100644 --- a/Layers/FullAtomModel/main.cpp +++ b/Layers/FullAtomModel/main.cpp @@ -16,6 +16,7 @@ PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) { m.def("getSeqNumAtoms", &getSeqNumAtoms, "Get number of atoms in a sequence"); m.def("PDB2CoordsUnordered", &PDB2CoordsUnordered, "Convert PDB to coordinates in the PDB order"); + m.def("PDB2CoordsOrdered", &PDB2CoordsOrdered, "Convert PDB to coordinates in the Angles2Coords order"); m.def("Coords2TypedCoords_forward", &Coords2TypedCoords_forward, "Convert coordinates to atom types"); m.def("Coords2TypedCoords_backward", &Coords2TypedCoords_backward, "Backward of Coords2TypedCoords"); diff --git a/Layers/Physics/AtomNames2Params/atomnames2params_interface.cpp b/Layers/Physics/AtomNames2Params/atomnames2params_interface.cpp new file mode 100644 index 0000000..c68c6aa --- /dev/null +++ b/Layers/Physics/AtomNames2Params/atomnames2params_interface.cpp @@ -0,0 +1,80 @@ +#include +#include "atomnames2params_interface.h" +#include "nUtil.h" +#include +#include + + +std::map AtomNames2Params_forward( torch::Tensor resnames, torch::Tensor atomnames, torch::Tensor num_atoms, + torch::Tensor types, torch::Tensor params, torch::Tensor assigned_params){ + CHECK_CPU_INPUT_TYPE(resnames, torch::kByte); + CHECK_CPU_INPUT_TYPE(atomnames, torch::kByte); + CHECK_CPU_INPUT_TYPE(types, torch::kByte); + CHECK_CPU_INPUT_TYPE(params, torch::kDouble); + CHECK_CPU_INPUT_TYPE(num_atoms, torch::kInt); + + int num_types = types.size(0); + std::map type_dict; + std::map assign_dict; + + auto param_acc = params.accessor(); + for(int i=0; i(); + + #pragma omp parallel for shared(type_dict, assign_dict) + for(int i=0; i &indexes){ + CHECK_CPU_INPUT_TYPE(gradOutput, torch::kDouble); + CHECK_CPU_INPUT_TYPE(gradInput, torch::kDouble); + + auto gradOutput_acc = gradOutput.accessor(); + + #pragma omg parallel for + for( auto element : indexes){ + int param_index = element.first; + assignment_indexes grad_indexes = element.second; + double grad_charge = 0.0; + double grad_radius = 0.0; + for(int i=0; i +#include +#include +#include + +typedef std::pair atom_type; +typedef std::pair atom_param; +typedef std::pair indexes_atom_param; +typedef std::vector> assignment_indexes; + +std::map AtomNames2Params_forward( torch::Tensor resnames, torch::Tensor atomnames, torch::Tensor num_atoms, + torch::Tensor types, torch::Tensor params, torch::Tensor assigned_types); + +void AtomNames2Params_backward( torch::Tensor gradOutput, torch::Tensor gradInput, + std::map &indexes); \ No newline at end of file diff --git a/Layers/Physics/Coords2Elec/coords2elec_interface.cpp b/Layers/Physics/Coords2Elec/coords2elec_interface.cpp new file mode 100644 index 0000000..b1526ec --- /dev/null +++ b/Layers/Physics/Coords2Elec/coords2elec_interface.cpp @@ -0,0 +1,83 @@ +#include "coords2elec_interface.h" +#include +#include "nUtil.h" +#include "KernelsElectrostatics.h" + +void Coords2Eps_forward(torch::Tensor coords, torch::Tensor assigned_params, torch::Tensor num_atoms, torch::Tensor eps, + float resolution, float ion_size, float wat_size, float asigma, int d){ + CHECK_GPU_INPUT_TYPE(coords, torch::kFloat); + CHECK_GPU_INPUT_TYPE(assigned_params, torch::kFloat); + CHECK_GPU_INPUT_TYPE(eps, torch::kFloat); + CHECK_GPU_INPUT_TYPE(num_atoms, torch::kInt); + + int batch_size = num_atoms.size(0); + #pragma omp parallel for + for(int i=0; i +#include "nUtil.h" +#include "KernelsStress.h" + +void Coords2Stress_forward( torch::Tensor coords, torch::Tensor vectors, torch::Tensor num_atoms, torch::Tensor volume, + float resolution){ + CHECK_GPU_INPUT_TYPE(coords, torch::kFloat); + CHECK_GPU_INPUT_TYPE(vectors, torch::kFloat); + CHECK_GPU_INPUT_TYPE(volume, torch::kFloat); + CHECK_GPU_INPUT_TYPE(num_atoms, torch::kInt); + + int batch_size = num_atoms.size(0); + #pragma omp parallel for + for(int i=0; i(), + single_vectors.data(), + num_atoms[i].item().toInt(), + single_volume.data(), + single_volume.size(1),//box_size + resolution); + } +} \ No newline at end of file diff --git a/Layers/Physics/Coords2Stress/coords2stress_interface.h b/Layers/Physics/Coords2Stress/coords2stress_interface.h new file mode 100644 index 0000000..46f5e25 --- /dev/null +++ b/Layers/Physics/Coords2Stress/coords2stress_interface.h @@ -0,0 +1,3 @@ +#include +void Coords2Stress_forward( torch::Tensor coords, torch::Tensor vectors, torch::Tensor num_atoms, torch::Tensor volume, + float resolution); \ No newline at end of file diff --git a/Layers/Physics/KernelsElectrostatics.cu b/Layers/Physics/KernelsElectrostatics.cu new file mode 100644 index 0000000..0decf14 --- /dev/null +++ b/Layers/Physics/KernelsElectrostatics.cu @@ -0,0 +1,231 @@ +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + + +__global__ void partialSumFaces( float* coords, float* assigned_params, int num_atoms, float* volume, + int box_size, float res, float ion_size, float wat_size, float asigma, uint d){ + + int vol_index = threadIdx.x; + float *single_volume = volume + vol_index * box_size*box_size*box_size; + + float cell_x, cell_y, cell_z; + float shift_cell_x=0.0, shift_cell_y=0.0, shift_cell_z=0.0; + float add_sigma = 0.0; + switch(vol_index){ + case 0: + shift_cell_x = 0.5; + add_sigma = wat_size; + break; + case 1: + shift_cell_y = 0.5; + add_sigma = wat_size; + break; + case 2: + shift_cell_z = 0.5; + add_sigma = wat_size; + break; + case 3: + add_sigma = ion_size; + break; + } + long cell_idx; + + float x, y, z; + int x_i, y_i, z_i; + float r2; + + + for(int idx = 0; idx=0 && i=0 && j=0 && k>>(coords, assigned_params, num_atoms, volume, box_size, res, ion_size, wat_size, asigma, d); +} + +void gpu_computeSumCells( float *coords, + float *assigned_params, + int num_atoms, + float *volume, + int box_size, + float res){ + sumCells<<<1, 1>>>(coords, assigned_params, num_atoms, volume, box_size, res); +} + +struct saxpy_functor +{ + const float mul; + saxpy_functor(float _mul):mul(_mul){} + + __host__ __device__ + float operator()(const float& x, const float& y) const + { + return mul * x + y; + } +}; + +struct fill_functor +{ + const size_t surf_size, box_size; + const float val; + fill_functor(size_t _surf_size, size_t _box_size, float _val) + :surf_size(_surf_size), box_size(_box_size), val(_val){} + + __host__ __device__ + float operator()(const size_t& idx, const float& old_val) const + { + int x = floor(idx/surf_size); + int y = floor((idx - x*surf_size)/box_size); + int z = idx - x*surf_size - y*box_size; + + if( (x==0||x==(box_size-1)) || (y==0||y==(box_size-1)) || (z==0||z==(box_size-1))){ + return val; + }else{ + return old_val; + } + } +}; + + +void gpu_computePhi( float *Q, float *Eps, float *Phi, size_t box_size, float res, float kappa02){ + size_t surf_size = box_size*box_size; + size_t vol_size = box_size*box_size*box_size; + size_t num_nonzero = 7*vol_size - 2*surf_size - 2*box_size - 2; + + cusp::dia_matrix A(vol_size, vol_size, num_nonzero, 7); + cusp::array1d phi(vol_size, 0.0); + cusp::array1d q(vol_size, 0.0); + A.diagonal_offsets[0] = -box_size*box_size; + A.diagonal_offsets[1] = -box_size; + A.diagonal_offsets[2] = -1; + A.diagonal_offsets[3] = 0; + A.diagonal_offsets[4] = 1; + A.diagonal_offsets[5] = box_size; + A.diagonal_offsets[6] = box_size*box_size; + + thrust::fill(A.values.values.begin(), A.values.values.end(), 0.0); + + thrust::device_ptr ei_begin(Eps); + thrust::device_ptr ej_begin(Eps + vol_size); + thrust::device_ptr ek_begin(Eps + 2*vol_size); + thrust::device_ptr lambda_begin(Eps + 3*vol_size); + thrust::device_ptr q_begin(Q); + thrust::device_ptr Phi_begin(Phi); + + //Lower diagonals + thrust::transform(ei_begin, ei_begin + vol_size - surf_size, A.values.column(0).begin() + surf_size, + A.values.column(0).begin() + surf_size, saxpy_functor(-res)); + thrust::transform(ej_begin, ej_begin + vol_size - box_size, A.values.column(1).begin() + box_size, + A.values.column(1).begin() + box_size, saxpy_functor(-res)); + thrust::transform(ek_begin, ek_begin + vol_size - 1, A.values.column(2).begin() + 1, + A.values.column(2).begin() + 1, saxpy_functor(-res)); + + //Upper diagonals + thrust::transform(ek_begin, ek_begin + vol_size - 1, A.values.column(4).begin(), + A.values.column(4).begin(), saxpy_functor(-res)); + thrust::transform(ej_begin, ej_begin + vol_size - box_size, A.values.column(5).begin(), + A.values.column(5).begin(), saxpy_functor(-res)); + thrust::transform(ei_begin, ei_begin + vol_size - surf_size, A.values.column(6).begin(), + A.values.column(6).begin(), saxpy_functor(-res)); + + //Diagonal + thrust::transform(ei_begin, ei_begin + vol_size, A.values.column(3).begin(), + A.values.column(3).begin(), saxpy_functor(res)); + thrust::transform(ej_begin, ej_begin + vol_size, A.values.column(3).begin(), + A.values.column(3).begin(), saxpy_functor(res)); + thrust::transform(ek_begin, ek_begin + vol_size, A.values.column(3).begin(), + A.values.column(3).begin(), saxpy_functor(res)); + + //diagonal shifted + thrust::transform(ei_begin, ei_begin + vol_size - surf_size, A.values.column(3).begin() + surf_size, + A.values.column(3).begin() + surf_size, saxpy_functor(res)); + thrust::transform(ej_begin, ej_begin + vol_size - box_size, A.values.column(3).begin() + box_size, + A.values.column(3).begin() + box_size, saxpy_functor(res)); + thrust::transform(ek_begin, ek_begin + vol_size - 1, A.values.column(3).begin() + 1, + A.values.column(3).begin() + 1, saxpy_functor(res)); + + //ionic term + thrust::transform(lambda_begin, lambda_begin + vol_size, A.values.column(3).begin(), + A.values.column(3).begin(), saxpy_functor(kappa02*res*res*res)); + + //boundary conditions (A_bound = 1) + thrust::counting_iterator first(0); + thrust::counting_iterator last(vol_size); + thrust::transform(first, last, A.values.column(3).begin(), A.values.column(3).begin(), fill_functor(surf_size, box_size, 1.0f)); + + thrust::transform(first, last, A.values.column(0).begin(), A.values.column(0).begin(), fill_functor(surf_size, box_size, 0.0f)); + thrust::transform(first, last, A.values.column(1).begin(), A.values.column(1).begin(), fill_functor(surf_size, box_size, 0.0f)); + thrust::transform(first, last, A.values.column(2).begin(), A.values.column(2).begin(), fill_functor(surf_size, box_size, 0.0f)); + + thrust::transform(first, last, A.values.column(4).begin(), A.values.column(4).begin(), fill_functor(surf_size, box_size, 0.0f)); + thrust::transform(first, last, A.values.column(5).begin(), A.values.column(5).begin(), fill_functor(surf_size, box_size, 0.0f)); + thrust::transform(first, last, A.values.column(6).begin(), A.values.column(6).begin(), fill_functor(surf_size, box_size, 0.0f)); + + + //charge + thrust::transform(q_begin, q_begin + vol_size, q.begin(), q.begin(), saxpy_functor(1.0)); + + cusp::monitor monitor(q, 1000, 1e-3, 0.0, false); + // monitor.set_verbose(); + cusp::precond::diagonal M(A); + cusp::krylov::cg(A, phi, q, monitor); + // monitor.print(); + + thrust::copy(phi.begin(), phi.end(), Phi); + + +} \ No newline at end of file diff --git a/Layers/Physics/KernelsElectrostatics.h b/Layers/Physics/KernelsElectrostatics.h new file mode 100644 index 0000000..f78500d --- /dev/null +++ b/Layers/Physics/KernelsElectrostatics.h @@ -0,0 +1,19 @@ + +void gpu_computePartialSumFaces( float *coords, + float *assigned_params, + int num_atoms, + float *volume, + int box_size, + float res, + float ion_size, float wat_size, float asigma, + uint d); + +void gpu_computeSumCells( float *coords, + float *assigned_params, + int num_atoms, + float *volume, + int box_size, + float res); + +void gpu_computePhi( float *Q, float *Eps, float *Phi, size_t box_size, + float res, float kappa02); \ No newline at end of file diff --git a/Layers/Physics/KernelsStress.cu b/Layers/Physics/KernelsStress.cu new file mode 100644 index 0000000..11261b7 --- /dev/null +++ b/Layers/Physics/KernelsStress.cu @@ -0,0 +1,53 @@ +#include +#include + + +__global__ void projectToGrid( float *coords, float *vectors, int num_atoms, float *volume, + int box_size, float res){ + int d=2; + int vol_index = threadIdx.x; + float *single_volume = volume + vol_index * box_size*box_size*box_size; + + float cell_x, cell_y, cell_z; + long cell_idx; + + float x, y, z; + int x_i, y_i, z_i; + float r2; + + + for(int idx = 0; idx=0 && i=0 && j=0 && k>>(coords, vectors, num_atoms, volume, box_size, res); +} \ No newline at end of file diff --git a/Layers/Physics/KernelsStress.h b/Layers/Physics/KernelsStress.h new file mode 100644 index 0000000..6e380f4 --- /dev/null +++ b/Layers/Physics/KernelsStress.h @@ -0,0 +1,6 @@ +void gpu_projectToGrid( float *coords, + float *vectors, + int num_atoms, + float *volume, + int box_size, + float res); \ No newline at end of file diff --git a/Layers/Physics/main.cpp b/Layers/Physics/main.cpp new file mode 100644 index 0000000..73948c0 --- /dev/null +++ b/Layers/Physics/main.cpp @@ -0,0 +1,18 @@ +#include +#include "coords2elec_interface.h" +#include "atomnames2params_interface.h" +#include "coords2stress_interface.h" +#include + +PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) { + m.def("AtomNames2Params_forward", &AtomNames2Params_forward, "Assigning parameters given atom names"); + m.def("AtomNames2Params_backward", &AtomNames2Params_backward, "Assigning parameters given atom names"); + + m.def("Coords2Eps_forward", &Coords2Eps_forward, "Projecting coordinates to the grid dielectric constant"); + m.def("Coords2Eps_backward", &Coords2Eps_backward, "Projecting coordinates to the grid dielectric constant"); + m.def("Coords2Q_forward", &Coords2Q_forward, "Projecting coordinates to the grid charge"); + m.def("QEps2Phi_forward", &QEps2Phi_forward, "Solving electrostatics"); + + m.def("Coords2Stress_forward", &Coords2Stress_forward, "Computing stress tensor from gamma"); + +} \ No newline at end of file diff --git a/Layers/RMSD/Coords2RMSD/coords2rmsd_interface.cpp b/Layers/RMSD/Coords2RMSD/coords2rmsd_interface.cpp index 4d138d7..6a7871b 100644 --- a/Layers/RMSD/Coords2RMSD/coords2rmsd_interface.cpp +++ b/Layers/RMSD/Coords2RMSD/coords2rmsd_interface.cpp @@ -124,7 +124,6 @@ void Coords2RMSD_forward( torch::Tensor centered_coords_src, if(centered_coords_src.ndimension() != 2){ ERROR("Incorrect input ndim"); } - int batch_size = num_atoms.size(0); for(int i=0; i -#include +#include "Coords2RMSD/coords2rmsd_interface.h" PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) { m.def("Coords2RMSDGPU_forward", &Coords2RMSDGPU_forward, "Coords2RMSD forward on GPU"); diff --git a/Layers/ReducedModel/Angles2Backbone/angles2backbone_interface.cpp b/Layers/ReducedModel/Angles2Backbone/angles2backbone_interface.cpp index fe19085..c094876 100644 --- a/Layers/ReducedModel/Angles2Backbone/angles2backbone_interface.cpp +++ b/Layers/ReducedModel/Angles2Backbone/angles2backbone_interface.cpp @@ -5,11 +5,13 @@ #include int Angles2BackboneGPU_forward(torch::Tensor input_angles, + torch::Tensor param, torch::Tensor output_coords, torch::Tensor angles_length, torch::Tensor A ){ CHECK_GPU_INPUT(input_angles); + CHECK_GPU_INPUT(param); CHECK_GPU_INPUT(output_coords); CHECK_GPU_INPUT(A); CHECK_GPU_INPUT_TYPE(angles_length, torch::kInt); @@ -18,22 +20,25 @@ int Angles2BackboneGPU_forward(torch::Tensor input_angles, } gpu_computeCoordinatesBackbone( input_angles.data(), + param.data(), output_coords.data(), A.data(), angles_length.data(), input_angles.size(0), input_angles.size(2)); } -int Angles2BackboneGPU_backward( torch::Tensor gradInput, - torch::Tensor gradOutput, - torch::Tensor input_angles, - torch::Tensor angles_length, - torch::Tensor A, - torch::Tensor dr_dangle +int Angles2BackboneGPUAngles_backward( torch::Tensor gradInput, + torch::Tensor gradOutput, + torch::Tensor input_angles, + torch::Tensor param, + torch::Tensor angles_length, + torch::Tensor A, + torch::Tensor dr_dangle ){ CHECK_GPU_INPUT(gradInput); CHECK_GPU_INPUT(gradOutput); CHECK_GPU_INPUT(input_angles); + CHECK_GPU_INPUT(param); CHECK_GPU_INPUT(dr_dangle); CHECK_GPU_INPUT(A); CHECK_GPU_INPUT_TYPE(angles_length, torch::kInt); @@ -41,6 +46,7 @@ int Angles2BackboneGPU_backward( torch::Tensor gradInput, ERROR("Incorrect input ndim"); } gpu_computeDerivativesBackbone( input_angles.data(), + param.data(), dr_dangle.data(), A.data(), angles_length.data(), @@ -55,14 +61,51 @@ int Angles2BackboneGPU_backward( torch::Tensor gradInput, input_angles.size(2)); } +int Angles2BackboneGPUParam_backward(torch::Tensor gradParam, + torch::Tensor gradOutput, + torch::Tensor input_angles, + torch::Tensor param, + torch::Tensor angles_length, + torch::Tensor A, + torch::Tensor dr_dparam + ){ + CHECK_GPU_INPUT(gradParam); + CHECK_GPU_INPUT(gradOutput); + CHECK_GPU_INPUT(input_angles); + CHECK_GPU_INPUT(param); + CHECK_GPU_INPUT(dr_dparam); + CHECK_GPU_INPUT(A); + CHECK_GPU_INPUT_TYPE(angles_length, torch::kInt); + if(gradOutput.ndimension() != 2){ + ERROR("Incorrect input ndim"); + } + + gpu_computeDerivativesParam( input_angles.data(), + param.data(), + dr_dparam.data(), + A.data(), + angles_length.data(), + input_angles.size(0), + input_angles.size(2)); + + gpu_backwardFromCoordsParam( gradParam.data(), + gradOutput.data(), + dr_dparam.data(), + angles_length.data(), + input_angles.size(0), + input_angles.size(2)); +} -int Angles2BackboneCPU_forward(torch::Tensor input_angles, - torch::Tensor output_coords, - torch::Tensor angles_length, - torch::Tensor A - ){ + +int Angles2BackboneCPU_forward( torch::Tensor input_angles, + torch::Tensor param, + torch::Tensor output_coords, + torch::Tensor angles_length, + torch::Tensor A + ){ CHECK_CPU_INPUT(input_angles); + CHECK_CPU_INPUT(param); CHECK_CPU_INPUT(output_coords); CHECK_CPU_INPUT(A); CHECK_CPU_INPUT_TYPE(angles_length, torch::kInt); @@ -71,17 +114,19 @@ int Angles2BackboneCPU_forward(torch::Tensor input_angles, } cpu_computeCoordinatesBackbone( input_angles.data(), - output_coords.data(), - A.data(), - angles_length.data(), - input_angles.size(0), - input_angles.size(2)); + param.data(), + output_coords.data(), + A.data(), + angles_length.data(), + input_angles.size(0), + input_angles.size(2)); } -int Angles2BackboneCPU_backward(torch::Tensor gradInput, +int Angles2BackboneCPUAngles_backward(torch::Tensor gradInput, torch::Tensor gradOutput, torch::Tensor input_angles, + torch::Tensor param, torch::Tensor angles_length, torch::Tensor A, torch::Tensor dr_dangle @@ -89,6 +134,7 @@ int Angles2BackboneCPU_backward(torch::Tensor gradInput, CHECK_CPU_INPUT(gradInput); CHECK_CPU_INPUT(gradOutput); CHECK_CPU_INPUT(input_angles); + CHECK_CPU_INPUT(param); CHECK_CPU_INPUT(dr_dangle); CHECK_CPU_INPUT(A); CHECK_CPU_INPUT_TYPE(angles_length, torch::kInt); @@ -97,6 +143,7 @@ int Angles2BackboneCPU_backward(torch::Tensor gradInput, } cpu_computeDerivativesBackbone( input_angles.data(), + param.data(), dr_dangle.data(), A.data(), angles_length.data(), @@ -109,4 +156,40 @@ int Angles2BackboneCPU_backward(torch::Tensor gradInput, angles_length.data(), input_angles.size(0), input_angles.size(2)); +} + +int Angles2BackboneCPUParam_backward(torch::Tensor gradParam, + torch::Tensor gradOutput, + torch::Tensor input_angles, + torch::Tensor param, + torch::Tensor angles_length, + torch::Tensor A, + torch::Tensor dr_dparam + ){ + CHECK_CPU_INPUT(gradParam); + CHECK_CPU_INPUT(gradOutput); + CHECK_CPU_INPUT(input_angles); + CHECK_CPU_INPUT(param); + CHECK_CPU_INPUT(dr_dparam); + CHECK_CPU_INPUT(A); + CHECK_CPU_INPUT_TYPE(angles_length, torch::kInt); + if(gradOutput.ndimension() != 2){ + ERROR("Incorrect input ndim"); + } + + + cpu_computeDerivativesParam( input_angles.data(), + param.data(), + dr_dparam.data(), + A.data(), + angles_length.data(), + input_angles.size(0), + input_angles.size(2)); + + cpu_backwardFromCoordsParam( gradParam.data(), + gradOutput.data(), + dr_dparam.data(), + angles_length.data(), + input_angles.size(0), + input_angles.size(2)); } \ No newline at end of file diff --git a/Layers/ReducedModel/Angles2Backbone/angles2backbone_interface.h b/Layers/ReducedModel/Angles2Backbone/angles2backbone_interface.h index 1e4b27b..dfe9b47 100644 --- a/Layers/ReducedModel/Angles2Backbone/angles2backbone_interface.h +++ b/Layers/ReducedModel/Angles2Backbone/angles2backbone_interface.h @@ -1,29 +1,48 @@ #include int Angles2BackboneGPU_forward(torch::Tensor input_angles, + torch::Tensor param, torch::Tensor output_coords, torch::Tensor angles_length, torch::Tensor A ); - -int Angles2BackboneGPU_backward( torch::Tensor gradInput, +int Angles2BackboneGPUAngles_backward( torch::Tensor gradInput, + torch::Tensor gradOutput, + torch::Tensor input_angles, + torch::Tensor param, + torch::Tensor angles_length, + torch::Tensor A, + torch::Tensor dr_dangle + ); +int Angles2BackboneGPUParam_backward(torch::Tensor gradParam, torch::Tensor gradOutput, torch::Tensor input_angles, + torch::Tensor param, torch::Tensor angles_length, - torch::Tensor A, - torch::Tensor dr_dangle + torch::Tensor A, + torch::Tensor dr_dparam ); int Angles2BackboneCPU_forward(torch::Tensor input_angles, - torch::Tensor output_coords, - torch::Tensor angles_length, - torch::Tensor A + torch::Tensor param, + torch::Tensor output_coords, + torch::Tensor angles_length, + torch::Tensor A ); -int Angles2BackboneCPU_backward( torch::Tensor gradInput, +int Angles2BackboneCPUAngles_backward(torch::Tensor gradInput, torch::Tensor gradOutput, torch::Tensor input_angles, + torch::Tensor param, torch::Tensor angles_length, torch::Tensor A, torch::Tensor dr_dangle + ); +int Angles2BackboneCPUParam_backward(torch::Tensor gradParam, + torch::Tensor gradOutput, + torch::Tensor input_angles, + torch::Tensor param, + torch::Tensor angles_length, + torch::Tensor A, + torch::Tensor dr_dparam ); \ No newline at end of file diff --git a/Layers/ReducedModel/cBackboneProteinCPUKernels.cpp b/Layers/ReducedModel/cBackboneProteinCPUKernels.cpp index 834b3bf..e3bd2bf 100644 --- a/Layers/ReducedModel/cBackboneProteinCPUKernels.cpp +++ b/Layers/ReducedModel/cBackboneProteinCPUKernels.cpp @@ -3,22 +3,22 @@ #include #include -#define R_CA_C 1.525 -#define R_C_N 1.330 -#define R_N_CA 1.460 - -#define CA_C_N (M_PI - 2.1186) -#define C_N_CA (M_PI - 1.9391) -#define N_CA_C (M_PI - 2.061) - template -void cpu_computeCoordinatesBackbone( T *angles, +void cpu_computeCoordinatesBackbone( T *angles, + T *param, T *atoms, T *A, int *length, int batch_size, int angles_stride){ int atoms_stride = 3*angles_stride; + + T R_N_CA = param[0]; + T C_N_CA = param[1]; + T R_CA_C = param[2]; + T N_CA_C = param[3]; + T R_C_N = param[4]; + T CA_C_N = param[5]; for(int batch_idx=0; batch_idx zero; zero.setZero(); cVector3 grad(dR_dAngle); - if( (3*angle_idx+angle_k) > atom_idx){ + if( (3*angle_idx+angle_k+1) > atom_idx){ grad.setZero(); }else{ @@ -88,10 +95,12 @@ Computes derivative of atom "atom_idx" coordinates with respect to {phi, psi, om grad = (Al * B * Ar_inv * Aj) * zero; } + } template -void cpu_computeDerivativesBackbone( T *angles, +void cpu_computeDerivativesBackbone( T *angles, + T *param, T *dR_dangle, T *A, int *length, @@ -106,10 +115,11 @@ void cpu_computeDerivativesBackbone( T *angles, for(int atom_idx=0; atom_idx( angles + (3*batch_idx+angle_k)*angles_stride, dR_dangle + (3*batch_idx+angle_k) * (atoms_stride*angles_stride*3) + angle_idx*(atoms_stride*3) + atom_idx*3, - d_A, angle_k, angle_idx, atom_idx); + d_A, angle_k, angle_idx, atom_idx, param); } } } @@ -144,14 +154,127 @@ void cpu_backwardFromCoordsBackbone( T *gradInput, } }; -template void cpu_computeCoordinatesBackbone( float*, float*, float*, int*, int, int); -template void cpu_computeCoordinatesBackbone( double*, double*, double*, int*, int, int); -template void device_singleAngleAtom( float*, float*, float*, int, int, int); -template void device_singleAngleAtom( double*, double*, double*, int, int, int); -template void cpu_computeDerivativesBackbone( float*, float*, float*, int*, int, int); -template void cpu_computeDerivativesBackbone( double*, double*, double*, int*, int, int); +template +void device_singleParamAtom(T *d_angle, //pointer to the angle stride + T *dR_dParamR, //pointer to the gradient + T *dR_dParamPsi, //pointer to the gradient + T *d_A, //pointer to the atom transformation matrix batch + int angle_k, //angle index (phi:0, psi:1, omega:2) + int angle_idx, //angle index + int atom_idx, //atom index + T *param +){ +/* +Computes derivative of atom "atom_idx" coordinates with respect to parameter in the matrix with the index "angle_idx". +*/ + T R_N_CA = param[0]; + T C_N_CA = param[1]; + T R_CA_C = param[2]; + T N_CA_C = param[3]; + T R_C_N = param[4]; + T CA_C_N = param[5]; + + T bond_angles[] = {C_N_CA, N_CA_C, CA_C_N}; + T bond_lengths[] = {R_N_CA, R_CA_C, R_C_N}; + cVector3 zero; zero.setZero(); + cVector3 gradR(dR_dParamR); + cVector3 gradPsi(dR_dParamPsi); + if( (3*angle_idx+angle_k+1) > atom_idx){ + gradR.setZero(); + gradPsi.setZero(); + }else{ + // std::cout<<"a"< Br, Bpsi, Ar_inv; + cMatrix44 Al(d_A + (3*angle_idx + angle_k)*16), Ar(d_A + (3*angle_idx + angle_k + 1)*16); + // std::cout<<"b"< Aj(d_A + atom_idx*16); + // std::cout<<"c"< +void cpu_computeDerivativesParam(T *angles, + T *param, + T *dR_dparam, + T *A, + int *length, + int batch_size, + int angles_stride){ + + int atoms_stride = 3*angles_stride; + for(int batch_idx=0; batch_idx( angles + (3*batch_idx+angle_k)*angles_stride, + dR_dparam + (6*batch_idx+2*angle_k)*(atoms_stride*angles_stride*3) + angle_idx*(atoms_stride*3) + atom_idx*3, + dR_dparam + (6*batch_idx+2*angle_k+1)*(atoms_stride*angles_stride*3) + angle_idx*(atoms_stride*3) + atom_idx*3, + d_A, angle_k, angle_idx, atom_idx, param); + } + } + } + } +}; + +template +void cpu_backwardFromCoordsParam( T *gradParam, + T *gradOutput, + T *dR_dparam, + int *length, + int batch_size, + int angles_stride){ + + int atoms_stride = 3*angles_stride; + for(int param_k=0; param_k<6; param_k++){ + for(int batch_idx=0; batch_idx dr(gradOutput + batch_idx*atoms_stride*3 + 3*atom_idx); + cVector3 dr_dangle(dR_dParam + atom_idx*3); + gradParam[param_k] += dr | dr_dangle; + // std::cout<( float*, float*, float*, float*, int*, int, int); +template void cpu_computeCoordinatesBackbone( double*, double*, double*, double*, int*, int, int); + +template void device_singleAngleAtom( float*, float*, float*, int, int, int, float*); +template void device_singleAngleAtom( double*, double*, double*, int, int, int, double*); + +template void cpu_computeDerivativesBackbone( float*, float*, float*, float*, int*, int, int); +template void cpu_computeDerivativesBackbone( double*, double*, double*, double*, int*, int, int); template void cpu_backwardFromCoordsBackbone( float*, float*, float*, int*, int, int); -template void cpu_backwardFromCoordsBackbone( double*, double*, double*, int*, int, int); \ No newline at end of file +template void cpu_backwardFromCoordsBackbone( double*, double*, double*, int*, int, int); + +template void cpu_computeDerivativesParam( float*, float*, float*, float*, int*, int, int); +template void cpu_computeDerivativesParam( double*, double*, double*, double*, int*, int, int); + +template void cpu_backwardFromCoordsParam( float*, float*, float*, int*, int, int); +template void cpu_backwardFromCoordsParam( double*, double*, double*, int*, int, int); \ No newline at end of file diff --git a/Layers/ReducedModel/cBackboneProteinCPUKernels.hpp b/Layers/ReducedModel/cBackboneProteinCPUKernels.hpp index a782105..8ef9720 100644 --- a/Layers/ReducedModel/cBackboneProteinCPUKernels.hpp +++ b/Layers/ReducedModel/cBackboneProteinCPUKernels.hpp @@ -1,5 +1,6 @@ template void cpu_computeCoordinatesBackbone( T *angles, + T *param, T *dr, T *dR_dangle, int *length, @@ -7,6 +8,7 @@ void cpu_computeCoordinatesBackbone( T *angles, int angles_stride); template void cpu_computeDerivativesBackbone( T *angles, + T *param, T *dR_dangle, T *A, int *length, @@ -18,4 +20,21 @@ void cpu_backwardFromCoordsBackbone( T *gradInput, T *dR_dangle, int *length, int batch_size, - int angles_stride); \ No newline at end of file + int angles_stride); + +template +void cpu_computeDerivativesParam( T *angles, + T *param, + T *dR_dparam, + T *A, + int *length, + int batch_size, + int angles_stride); + +template +void cpu_backwardFromCoordsParam( T *gradParam, + T *gradOutput, + T *dR_dparam, + int *length, + int batch_size, + int angles_stride); \ No newline at end of file diff --git a/Layers/ReducedModel/cBackboneProteinCUDAKernels.cu b/Layers/ReducedModel/cBackboneProteinCUDAKernels.cu index 3c4b2b2..15800e5 100644 --- a/Layers/ReducedModel/cBackboneProteinCUDAKernels.cu +++ b/Layers/ReducedModel/cBackboneProteinCUDAKernels.cu @@ -4,11 +4,18 @@ #define WARP_SIZE 32 template -__global__ void computeCoordinatesBackbone( T *angles, T *atoms, T *A, int *length, int angles_stride){ +__global__ void computeCoordinatesBackbone( T *angles, T *param, T *atoms, T *A, int *length, int angles_stride){ uint batch_idx = threadIdx.x; int atoms_stride = 3*angles_stride; int num_angles = length[batch_idx]; int num_atoms = 3*length[batch_idx]; + + T R_N_CA = param[0]; + T C_N_CA = param[1]; + T R_CA_C = param[2]; + T N_CA_C = param[3]; + T R_C_N = param[4]; + T CA_C_N = param[5]; // printf("angles_stride=%d, batch_idx=%d, num_atoms=%d, num_angles=%d\n",angles_stride, batch_idx,num_atoms,num_angles); T *d_atoms = atoms + batch_idx*(atoms_stride)*3; @@ -39,64 +46,31 @@ __global__ void computeCoordinatesBackbone( T *angles, T *atoms, T *A, int *leng // printf("batch_idx=%d, i=%d, x=%f, y=%f, z=%f\n", batch_idx, i, *(d_atoms + 3*i), *(d_atoms + 3*i+1), *(d_atoms + 3*i+2)); } } -template -__global__ void computeGradientsOptimizedBackboneOmega(T *angles, T *dR_dangle, T *A, int *length, int angles_stride){ - - uint batch_size = blockDim.x; - uint batch_idx = blockIdx.x; - uint atom_i_idx = blockIdx.y; - uint local_angle_k_idx = threadIdx.x; - uint angle_k_idx = blockIdx.z*WARP_SIZE + local_angle_k_idx; - - int num_angles = length[batch_idx]; - int num_atoms = 3*num_angles; - int atoms_stride = 3*angles_stride; - - T *d_A = A + batch_idx * atoms_stride * 16; - T *d_A_warp = d_A + 3*blockIdx.z*WARP_SIZE*16; - __shared__ T s_A [WARP_SIZE*3*16]; - for(int i=local_angle_k_idx; i=angles_stride)return; - T *d_omega = angles + (3*batch_idx+2)*angles_stride; - T *dR_dOmega = dR_dangle + (3*batch_idx+2) * (atoms_stride*angles_stride*3) + atom_i_idx*angles_stride*3 + angle_k_idx*3; - T tmp1[16], tmp2[16], tmp3[16]; - - //dA_i / dpsi_k - if( (3*angle_k_idx) > atom_i_idx || (3*angle_k_idx == 0)){ - setVec3(dR_dOmega, 0, 0, 0); - }else{ - getRotationMatrixDihedralDPsi(tmp2, d_omega[angle_k_idx], CA_C_N, R_C_N); - mat44Mul(s_A + (3*local_angle_k_idx-1)*16, tmp2, tmp1); - - invertMat44(tmp3, s_A + (3*local_angle_k_idx)*16); - mat44Mul(tmp3, d_A + atom_i_idx*16, tmp2); - - mat44Mul(tmp1, tmp2, tmp3); - mat44Zero3Mul(tmp3, dR_dOmega); - } - -} template __device__ void device_singleAngleAtom( T *d_angle, //pointer to the angle stride - T *dR_dAngle, //pointer to the gradient + T *dR_dAngle, //pointer to the gradient of angles T *d_A, //pointer to the atom transformation matrix batch T *d_A_warp, //shared pointer to the atom transformation matrix int angle_k, //angle index (phi:0, psi:1, omega:2) int angle_idx, //angle index int angle_widx, //angle index in the warp - int atom_idx //atom index + int atom_idx, //atom index + T *param ){ + T R_N_CA = param[0]; + T C_N_CA = param[1]; + T R_CA_C = param[2]; + T N_CA_C = param[3]; + T R_C_N = param[4]; + T CA_C_N = param[5]; + T tmp1[16], tmp2[16], tmp3[16]; T B[16], Ar_inv[16]; T bond_angles[] = {C_N_CA, N_CA_C, CA_C_N}; T bond_lengths[] = {R_N_CA, R_CA_C, R_C_N}; - if( (3*angle_idx+angle_k) > atom_idx){ + if( (3*angle_idx+angle_k+1) > atom_idx){ setVec3(dR_dAngle, 0, 0, 0); }else{ getRotationMatrixDihedralDPsi(B, d_angle[angle_idx], bond_angles[angle_k], bond_lengths[angle_k]); @@ -110,7 +84,7 @@ __device__ void device_singleAngleAtom( T *d_angle, //pointer to the angle strid } } template -__global__ void computeGradientsOptimizedBackbone( T *angles, T *dR_dangle, T *A, int *length, int angles_stride){ +__global__ void computeGradientsOptimizedBackbone( T *angles, T *param, T *dR_dangle, T *A, int *length, int angles_stride){ uint batch_idx = blockIdx.x; uint atom_idx = blockIdx.y; uint angle_widx = threadIdx.x; @@ -139,7 +113,7 @@ __global__ void computeGradientsOptimizedBackbone( T *angles, T *dR_dangle, T *A dR_dangle + (3*batch_idx + angle_k) * (atoms_stride*angles_stride*3) + angle_idx*atoms_stride*3 + atom_idx*3, d_A, s_A, - angle_k, angle_idx, angle_widx, atom_idx); + angle_k, angle_idx, angle_widx, atom_idx, param); } } @@ -162,14 +136,14 @@ __global__ void backwardFromCoordinatesBackbone(T *angles, T *dr, T *dR_dangle, } } template -void gpu_computeCoordinatesBackbone(T *angles, T *atoms, T *A, int *length, int batch_size, int angles_stride){ - computeCoordinatesBackbone<<<1, batch_size>>>(angles, atoms, A, length, angles_stride); +void gpu_computeCoordinatesBackbone(T *angles, T *param, T *atoms, T *A, int *length, int batch_size, int angles_stride){ + computeCoordinatesBackbone<<<1, batch_size>>>(angles, param, atoms, A, length, angles_stride); } template -void gpu_computeDerivativesBackbone(T *angles, T *dR_dangle, T *A, int *length, int batch_size, int angles_stride){ +void gpu_computeDerivativesBackbone(T *angles, T *param, T *dR_dangle, T *A, int *length, int batch_size, int angles_stride){ dim3 batch_angles_dim_special(batch_size, 3*angles_stride, angles_stride/WARP_SIZE + 1); - computeGradientsOptimizedBackbone<<>>(angles, dR_dangle, A, length, angles_stride); + computeGradientsOptimizedBackbone<<>>(angles, param, dR_dangle, A, length, angles_stride); } template void gpu_backwardFromCoordsBackbone(T *angles, T *dr, T *dR_dangle, int *length, int batch_size, int angles_stride){ @@ -177,10 +151,133 @@ void gpu_backwardFromCoordsBackbone(T *angles, T *dr, T *dR_dangle, int *length, backwardFromCoordinatesBackbone<<>>(angles, dr, dR_dangle, length, angles_stride); } -template void gpu_computeCoordinatesBackbone(float*, float*, float*, int*, int, int); -template void gpu_computeDerivativesBackbone(float*, float*, float*, int*, int, int); + + +template +__device__ void device_singleParamAtom( T *d_angle, //pointer to the angle stride + T *dR_dParamR, //pointer to the gradient of parameters + T *dR_dParamPsi, //pointer to the gradient of parameters + T *d_A, //pointer to the atom transformation matrix batch + T *d_A_warp, //shared pointer to the atom transformation matrix + int angle_k, //angle index (phi:0, psi:1, omega:2) + int angle_idx, //angle index + int angle_widx, //angle index in the warp + int atom_idx, //atom index + T *param +){ + T R_N_CA = param[0]; + T C_N_CA = param[1]; + T R_CA_C = param[2]; + T N_CA_C = param[3]; + T R_C_N = param[4]; + T CA_C_N = param[5]; + + T tmp1r[16], tmp1psi[16], tmp2[16], tmp3[16]; + T Br[16], Bpsi[16], Ar_inv[16]; + T bond_angles[] = {C_N_CA, N_CA_C, CA_C_N}; + T bond_lengths[] = {R_N_CA, R_CA_C, R_C_N}; + + if( (3*angle_idx+angle_k+1) > atom_idx){ + setVec3(dR_dParamR, 0, 0, 0); + setVec3(dR_dParamPsi, 0, 0, 0); + }else{ + getRotationMatrixDihedralDr(Br, d_angle[angle_idx], bond_angles[angle_k], bond_lengths[angle_k]); + getRotationMatrixDihedralDkappa(Bpsi, d_angle[angle_idx], bond_angles[angle_k], bond_lengths[angle_k]); + mat44Mul(d_A_warp + (3*angle_widx + angle_k)*16, Br, tmp1r); + mat44Mul(d_A_warp + (3*angle_widx + angle_k)*16, Bpsi, tmp1psi); + + invertMat44(Ar_inv, d_A_warp + (3*angle_widx + angle_k + 1)*16); + mat44Mul(Ar_inv, d_A + atom_idx*16, tmp2); + + mat44Mul(tmp1r, tmp2, tmp3); + mat44Zero3Mul(tmp3, dR_dParamR); + mat44Mul(tmp1psi, tmp2, tmp3); + mat44Zero3Mul(tmp3, dR_dParamPsi); + // printf("%d %d %d: (%f %f %f), (%f %f %f)\n", angle_idx, angle_k, atom_idx, *(dR_dParamR), *(dR_dParamR+1), *(dR_dParamR+2), + // *(dR_dParamPsi), *(dR_dParamPsi+1), *(dR_dParamPsi+2)); + + // T *A = d_A_warp + (3*angle_widx + angle_k)*16; + // T *A = d_A + atom_idx*16; + // printf("%d:\n%f %f %f\n%f %f %f\n%f %f %f\n",atom_idx*16, A[0], A[1], A[2], A[4], A[5], A[6], A[8], A[9], A[10]); + } +} +template +__global__ void computeGradientsOptimizedParam( T *angles, T *param, T *dR_dparam, T *A, int *length, int angles_stride){ + uint batch_idx = blockIdx.x; + uint atom_idx = blockIdx.y; + uint angle_widx = threadIdx.x; + uint warp_idx = blockIdx.z; + uint angle_idx = warp_idx*WARP_SIZE + angle_widx; + + int num_angles = length[batch_idx]; + int num_atoms = 3*num_angles; + int atoms_stride = 3*angles_stride; + + T *d_A = A + batch_idx * atoms_stride * 16; + T *d_A_warp = d_A + 3*warp_idx*WARP_SIZE * 16; + __shared__ T s_A[(3*WARP_SIZE + 1)*16]; + + for(int i=threadIdx.x; i<(3*WARP_SIZE + 1)*16; i+=WARP_SIZE){ + if( ((3*warp_idx*WARP_SIZE + 1)*16 + i) < (num_atoms*16)){ + s_A[i] = d_A_warp[i]; + } + } + __syncthreads(); + + if(angle_idx>=angles_stride)return; + + for(int angle_k=0; angle_k<3; angle_k++){ + // printf("%d, %d, %d, %d: %d\n",batch_idx, atom_idx, angle_idx, angle_k, + // (6*batch_idx + 2*angle_k) * (atoms_stride*angles_stride*3) + angle_idx*atoms_stride*3 + atom_idx*3); + device_singleParamAtom( angles + (3*batch_idx + angle_k) * angles_stride, + dR_dparam + (6*batch_idx + 2*angle_k) * (atoms_stride*angles_stride*3) + angle_idx*atoms_stride*3 + atom_idx*3, + dR_dparam + (6*batch_idx + 2*angle_k+1) * (atoms_stride*angles_stride*3) + angle_idx*atoms_stride*3 + atom_idx*3, + d_A, + s_A, + angle_k, angle_idx, angle_widx, atom_idx, param); + } + +} +template +__global__ void backwardFromCoordinatesParam(T *gradParam, T *gradOutput, T *dR_dparam, int *length, int angles_stride, int batch_size){ + + int param_k = threadIdx.x; + int atoms_stride = 3*angles_stride; + for(int batch_idx=0; batch_idx(gradOutput + batch_idx*atoms_stride*3 + 3*atom_idx, dR_dParam + 3*atom_idx); + } + } + } +} +template +void gpu_computeDerivativesParam(T *angles, T *param, T *dR_dparam, T *A, int *length, int batch_size, int angles_stride){ + dim3 batch_angles_dim_special(batch_size, 3*angles_stride, angles_stride/WARP_SIZE + 1); + computeGradientsOptimizedParam<<>>(angles, param, dR_dparam, A, length, angles_stride); +} +template +void gpu_backwardFromCoordsParam(T *gradParam, T *dr, T *dR_dparam, int *length, int batch_size, int angles_stride){ + backwardFromCoordinatesParam<<<1, 6>>>(gradParam, dr, dR_dparam, length, angles_stride, batch_size); +} + +template void gpu_computeCoordinatesBackbone(float*, float*, float*, float*, int*, int, int); +template void gpu_computeDerivativesBackbone(float*, float*, float*, float*, int*, int, int); template void gpu_backwardFromCoordsBackbone(float*, float*, float*, int*, int, int); -template void gpu_computeCoordinatesBackbone(double*, double*, double*, int*, int, int); -template void gpu_computeDerivativesBackbone(double*, double*, double*, int*, int, int); +template void gpu_computeDerivativesParam(float*, float*, float*, float*, int*, int, int); +template void gpu_backwardFromCoordsParam(float*, float*, float*, int*, int, int); + + +template void gpu_computeCoordinatesBackbone(double*, double*, double*, double*, int*, int, int); +template void gpu_computeDerivativesBackbone(double*, double*, double*, double*, int*, int, int); template void gpu_backwardFromCoordsBackbone(double*, double*, double*, int*, int, int); + +template void gpu_computeDerivativesParam(double*, double*, double*, double*, int*, int, int); +template void gpu_backwardFromCoordsParam(double*, double*, double*, int*, int, int); diff --git a/Layers/ReducedModel/cBackboneProteinCUDAKernels.h b/Layers/ReducedModel/cBackboneProteinCUDAKernels.h index c30dc5b..f753de1 100644 --- a/Layers/ReducedModel/cBackboneProteinCUDAKernels.h +++ b/Layers/ReducedModel/cBackboneProteinCUDAKernels.h @@ -1,23 +1,41 @@ -#define REAL float template void gpu_computeCoordinatesBackbone(T *angles, - T *dr, - T *dR_dangle, - int *length, - int batch_size, - int angles_stride); + T *param, + T *atoms, + T *A, + int *length, + int batch_size, + int angles_stride); +template +void gpu_computeDerivativesBackbone(T *angles, + T *param, + T *dR_dangle, + T *A, + int *length, + int batch_size, + int angles_stride); template -void gpu_computeDerivativesBackbone(T *angles, - T *atoms, - T *A, - int *length, - int batch_size, +void gpu_backwardFromCoordsBackbone(T *angles, + T *dr, + T *dR_dangle, + int *length, + int batch_size, int angles_stride); + +template +void gpu_computeDerivativesParam( T *angles, + T *param, + T *dR_dparam, + T *A, + int *length, + int batch_size, + int angles_stride); + template -void gpu_backwardFromCoordsBackbone(T *angles, - T *dR_dangle, - T *A, - int *length, - int batch_size, +void gpu_backwardFromCoordsParam( T *gradParam, + T *dr, + T *dR_dparam, + int *length, + int batch_size, int angles_stride); \ No newline at end of file diff --git a/Layers/ReducedModel/main.cpp b/Layers/ReducedModel/main.cpp index e78d8e5..717e957 100644 --- a/Layers/ReducedModel/main.cpp +++ b/Layers/ReducedModel/main.cpp @@ -3,7 +3,9 @@ PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) { m.def("Angles2BackboneCPU_forward", &Angles2BackboneCPU_forward, "Angles2Backbone cpu forward"); - m.def("Angles2BackboneCPU_backward", &Angles2BackboneCPU_backward, "Angles2Backbone cpu backward"); + m.def("Angles2BackboneCPUAngles_backward", &Angles2BackboneCPUAngles_backward, "Angles2Backbone cpu backward for angles"); + m.def("Angles2BackboneCPUParam_backward", &Angles2BackboneCPUParam_backward, "Angles2Backbone cpu backward for parameters"); m.def("Angles2BackboneGPU_forward", &Angles2BackboneGPU_forward, "Angles2Backbone gpu forward"); - m.def("Angles2BackboneGPU_backward", &Angles2BackboneGPU_backward, "Angles2Backbone gpu backward"); + m.def("Angles2BackboneGPUAngles_backward", &Angles2BackboneGPUAngles_backward, "Angles2Backbone gpu backward for angles"); + m.def("Angles2BackboneGPUParam_backward", &Angles2BackboneGPUParam_backward, "Angles2Backbone gpu backward for parameters"); } \ No newline at end of file diff --git a/Layers/Volume/HashKernel.cu b/Layers/Volume/HashKernel.cu new file mode 100644 index 0000000..c6f7733 --- /dev/null +++ b/Layers/Volume/HashKernel.cu @@ -0,0 +1,267 @@ +#include +#include +#include +namespace cg = cooperative_groups; + + +#include "thrust/device_ptr.h" +#include "thrust/for_each.h" +#include "thrust/iterator/zip_iterator.h" +#include "thrust/sort.h" + +#include "HashKernel.h" +#include "cub/cub.cuh" + +#define HASH_EMPTY 2147483647 + +//Round a / b to nearest higher integer value +uint iDivUp(uint a, uint b) +{ + return (a % b != 0) ? (a / b + 1) : (a / b); +} + +// compute grid and thread block size for a given number of elements +void computeGridSize(uint n, uint blockSize, uint &numBlocks, uint &numThreads) +{ + numThreads = min(blockSize, n); + numBlocks = iDivUp(n, numThreads); +} + +template +__device__ int3 computeGridPos(T p_x, T p_y, T p_z, float res){ + int3 gridPos; + gridPos.x = floor(p_x/res); + gridPos.y = floor(p_y/res); + gridPos.z = floor(p_z/res); + return gridPos; +} + +__device__ long computeGridHash(int3 gridPos, int spatial_dim){ + return gridPos.z + gridPos.y * spatial_dim + gridPos.x * spatial_dim * spatial_dim; +} + +template +__global__ void computeHash(long *gridParticleHash, long *gridParticleIndex, T *coords, + int num_atoms, int spatial_dim, float res, int d){ + uint atom_idx = blockIdx.x*blockDim.x + threadIdx.x; + uint num_neighbours = (2*d+1)*(2*d+1)*(2*d+1); + if (atom_idx>=num_atoms) return; + + int3 gridPos = computeGridPos(coords[3*atom_idx+0], coords[3*atom_idx+1], coords[3*atom_idx+2], res); + uint neighbour_idx = 0; + for(int i=gridPos.x-d; i<=gridPos.x+d; i++){ + for(int j=gridPos.y-d; j<=gridPos.y+d; j++){ + for(int k=gridPos.z-d; k<=gridPos.z+d; k++){ + if( ((i>=0 && i=0 && j=0 && k +__global__ void reorderData(long *cellStart, long *cellEnd, + T *sortedPos, + long *gridParticleHash, long *gridParticleIndex, + T *coords, int num_atoms){ + + cg::thread_block cta = cg::this_thread_block(); + extern __shared__ long sharedHash[]; + uint index = blockIdx.x * blockDim.x + threadIdx.x; + long hash; + + if(index < num_atoms){ + hash = gridParticleHash[index]; + sharedHash[threadIdx.x + 1] = hash; + if(threadIdx.x == 0 && index > 0){ + sharedHash[0] = gridParticleHash[index-1]; + } + + } + + cg::sync(cta); + + if(index < num_atoms){ + + if(sharedHash[threadIdx.x] == HASH_EMPTY){ + return; + } + + if(index == 0 || hash != sharedHash[threadIdx.x]){ + if(hash != HASH_EMPTY) + cellStart[hash] = index; + if (index > 0) + cellEnd[sharedHash[threadIdx.x]] = index; + } + if(index == num_atoms - 1 && hash != HASH_EMPTY){ + cellEnd[hash] = index + 1; + } + if(hash != HASH_EMPTY){ + long sortedIndex = gridParticleIndex[index]; + sortedPos[3*index+0] = coords[3*sortedIndex+0]; + sortedPos[3*index+1] = coords[3*sortedIndex+1]; + sortedPos[3*index+2] = coords[3*sortedIndex+2]; + } + } +} + + +template +__global__ void evalCell(T *sortedPos, T *volume, + long *cellStart, long *cellEnd, + int spatial_dim, float res){ + + uint i = (blockIdx.x * blockDim.x) + threadIdx.x; + uint j = (blockIdx.y * blockDim.y) + threadIdx.y; + uint k = (blockIdx.z * blockDim.z) + threadIdx.z; + + if( !((i>=0 && i=0 && j=0 && k +void gpu_computeCoords2Volume( T *coords, + int num_atoms, + T *volume, + int spatial_dim, + float res, + int d, + long *gridParticleHash, + long *gridParticleHash_buf, + long *gridParticleIndex, + long *gridParticleIndex_buf, + long *cellStart, + long *cellStop, + T *sortedPos, + int dev){ + cudaSetDevice(dev); + if(num_atoms == 0)return; + uint num_neighbours = (2*d+1)*(2*d+1)*(2*d+1); + + uint numThreads, numBlocks; + computeGridSize(num_atoms, 64, numBlocks, numThreads); + computeHash<<>>( gridParticleHash, gridParticleIndex, + coords, num_atoms, spatial_dim, res, d); + gpuErrchk( cudaPeekAtLastError() ); + + void *d_temp_storage = NULL; + size_t temp_storage_bytes = 0; + cub::DoubleBuffer b_gridParticleHash(gridParticleHash, gridParticleHash_buf); + cub::DoubleBuffer b_gridParticleIndex(gridParticleIndex, gridParticleIndex_buf); + cub::DeviceRadixSort::SortPairs(d_temp_storage, temp_storage_bytes, b_gridParticleHash, b_gridParticleIndex, num_neighbours*num_atoms); + cudaMalloc(&d_temp_storage, temp_storage_bytes); + cub::DeviceRadixSort::SortPairs(d_temp_storage, temp_storage_bytes, b_gridParticleHash, b_gridParticleIndex, num_neighbours*num_atoms); + cudaFree(d_temp_storage); + gpuErrchk( cudaPeekAtLastError() ); + + + computeGridSize(num_neighbours*num_atoms, 64, numBlocks, numThreads); + uint smemSize = sizeof(long)*(numThreads+1); + reorderData<<>>(cellStart, cellStop, + sortedPos, + b_gridParticleHash.Current(), b_gridParticleIndex.Current(), + coords, num_neighbours*num_atoms); + gpuErrchk( cudaPeekAtLastError() ); + + dim3 d3ThreadsPerBlock(4, 4, 4); + dim3 d3NumBlocks( spatial_dim/d3ThreadsPerBlock.x + 1, + spatial_dim/d3ThreadsPerBlock.y + 1, + spatial_dim/d3ThreadsPerBlock.z + 1); + evalCell<<>>(sortedPos, volume, cellStart, cellStop, spatial_dim, res); + + // cudaFree(gridParticleHash); + // cudaFree(gridParticleIndex); + // cudaFree(cellStart); + // cudaFree(cellEnd); + // cudaFree(sortedPos); +} + + + +template +__global__ void projectFromTensor( T* coords, T* grad, int num_atoms, + T *volume, int spatial_dim, float res, int d){ + uint atom_idx = (blockIdx.x * blockDim.x) + threadIdx.x; + if(atom_idx>num_atoms)return; + + T x = coords[3*atom_idx], y = coords[3*atom_idx + 1], z = coords[3*atom_idx + 2]; + int x_i = floor(x/res); + int y_i = floor(y/res); + int z_i = floor(z/res); + + T grad_x=0, grad_y=0, grad_z=0; + for(int i=x_i-d; i<=(x_i+d);i++){ + for(int j=y_i-d; j<=(y_i+d);j++){ + for(int k=z_i-d; k<=(z_i+d);k++){ + if( (i>=0 && i=0 && j=0 && k +void gpu_computeVolume2Coords( T *coords, + T* grad, + int num_atoms, + T *volume, + int spatial_dim, float res, int d){ + dim3 threadsPerBlock(64); + dim3 numBlocks( num_atoms/threadsPerBlock.x + 1); + + projectFromTensor<<>>(coords, grad, num_atoms, volume, spatial_dim, res, d); +} + + +template void gpu_computeVolume2Coords( float*, float*, int, float*, int, float, int); +template void gpu_computeVolume2Coords( double*, double*, int, double*, int, float, int); + +template void gpu_computeCoords2Volume(float*, int, float*, int, float, int, long*, long*, long*, long*, long*, long*, float*, int); +template void gpu_computeCoords2Volume(double*, int, double*, int, float, int, long*, long*, long*, long*, long*, long*, double*, int); diff --git a/Layers/Volume/HashKernel.h b/Layers/Volume/HashKernel.h new file mode 100644 index 0000000..d442f86 --- /dev/null +++ b/Layers/Volume/HashKernel.h @@ -0,0 +1,23 @@ +template +void gpu_computeCoords2Volume( T *coords, + int num_atoms, + T *volume, + int spatial_dim, + float res, int d, + long *gridParticleHash, + long *gridParticleHash_buf, + long *gridParticleIndex, + long *gridParticleIndex_buf, + long *cellStart, + long *cellStop, + T *sortedPos, + int dev); + +template +void gpu_computeVolume2Coords( T *coords, + T* grad, + int num_atoms, + T *volume, + int spatial_dim, + float res, + int d); \ No newline at end of file diff --git a/Layers/Volume/Kernels.cu b/Layers/Volume/Kernels.cu index eec3ae2..66febe5 100644 --- a/Layers/Volume/Kernels.cu +++ b/Layers/Volume/Kernels.cu @@ -1,99 +1,102 @@ #include -template -__global__ void projectToTensor(T* coords, int* num_atoms_of_type, int* offsets, T *volume, - int spatial_dim, float res){ /* -Input: - coords: coordinates in a flat array: - coords: {protein1, ... proteinN} - protein1: {atom_type1 .. atom_typeM} - atom_type: {x1,y1,z1 .. xL,yL,zL} - num_atoms_of_type: number of atoms in each atom_type - offsets: offset for coordinates for each atom_type volume -Output: - volume: density -*/ - int d = 2; - int type_index = threadIdx.x; - T *type_volume = volume + type_index * spatial_dim*spatial_dim*spatial_dim; - T *atoms_coords = coords + 3*offsets[type_index]; - int n_atoms = num_atoms_of_type[type_index]; - for(int atom_idx = 0; atom_idx<3*n_atoms; atom_idx+=3){ - T x = atoms_coords[atom_idx], - y = atoms_coords[atom_idx + 1], - z = atoms_coords[atom_idx + 2]; - int x_i = floor(x/res); +template +__global__ void projectToCell(T* coords, int num_atoms, T *volume, int spatial_dim, float res){ + uint d = 2; + uint i = (blockIdx.x * blockDim.x) + threadIdx.x; + uint j = (blockIdx.y * blockDim.y) + threadIdx.y; + uint k = (blockIdx.z * blockDim.z) + threadIdx.z; + + if( !((i>=0 && i=0 && j=0 && k=0 && i=0 && j=0 && kd) continue; + if(__sad(y_i, j, 0)>d) continue; + if(__sad(z_i, k, 0)>d) continue; + + T r2 = (x - i*res)*(x - i*res)+(y - j*res)*(y - j*res)+(z - k*res)*(z - k*res); + + result += exp(-r2/2.0); } + volume[idx] = result; } + template -__global__ void projectFromTensor(T* coords, T* grad, int* num_atoms_of_type, int* offsets, T *volume, - int spatial_dim, float res){ -/* -Input: - coords: coordinates in a flat array: - coords: {protein1, ... proteinN} - protein1: {atom_type1 .. atom_typeM} - atom_type: {x1,y1,z1 .. xL,yL,zL} - num_atoms_of_type: number of atoms in each atom_type - offsets: offset for coordinates for each atom_type volume - volume: gradient to be projected on atoms -Output: - grad: for each atom to store the gradient projection -*/ - int d = 2; - int type_index = threadIdx.x; - T *type_volume = volume + type_index * spatial_dim*spatial_dim*spatial_dim; - T *atoms_coords = coords + 3*offsets[type_index]; - T *grad_coords = grad + 3*offsets[type_index]; - int n_atoms = num_atoms_of_type[type_index]; - for(int atom_idx = 0; atom_idx<3*n_atoms; atom_idx+=3){ - T x = atoms_coords[atom_idx], - y = atoms_coords[atom_idx + 1], - z = atoms_coords[atom_idx + 2]; - // grad_coords[atom_idx] = 0.0; - // grad_coords[atom_idx+1] = 0.0; - // grad_coords[atom_idx+2] = 0.0; - int x_i = floor(x/res); - int y_i = floor(y/res); - int z_i = floor(z/res); - - for(int i=x_i-d; i<=(x_i+d);i++){ - for(int j=y_i-d; j<=(y_i+d);j++){ - for(int k=z_i-d; k<=(z_i+d);k++){ - if( (i>=0 && i=0 && j=0 && k<<>>(coords, num_atoms, volume, spatial_dim, res); + +} + + +template +__global__ void projectFromTensor( T* coords, T* grad, int num_atoms, + T *volume, int spatial_dim, float res){ + uint d = 2; + uint atom_idx = (blockIdx.x * blockDim.x) + threadIdx.x; + if(atom_idx>num_atoms)return; + + T x = coords[3*atom_idx], y = coords[3*atom_idx + 1], z = coords[3*atom_idx + 2]; + int x_i = floor(x/res); + int y_i = floor(y/res); + int z_i = floor(z/res); + + T grad_x=0, grad_y=0, grad_z=0; + for(int i=x_i-d; i<=(x_i+d);i++){ + for(int j=y_i-d; j<=(y_i+d);j++){ + for(int k=z_i-d; k<=(z_i+d);k++){ + if( (i>=0 && i=0 && j=0 && k +void gpu_computeVolume2Coords( T *coords, + T* grad, + int num_atoms, + T *volume, + int spatial_dim, + float res){ + dim3 threadsPerBlock(64); + dim3 numBlocks( num_atoms/threadsPerBlock.x + 1); -__global__ void selectFromTensor(float *features, + projectFromTensor<<>>(coords, grad, num_atoms, volume, spatial_dim, res); +} +*/ + +__global__ void coordSelect(float *features, float* volume, int spatial_dim, float *coords, int num_atoms, int max_num_atoms, float res){ @@ -122,47 +125,61 @@ Output: } } -template -void gpu_computeCoords2Volume( T *coords, - int *num_atoms_of_type, - int *offsets, - T *volume, - int spatial_dim, - int num_atom_types, +__global__ void coordSelectGrad(float *gradOutput, + float* gradInput, int spatial_dim, + float *coords, int num_atoms, int max_num_atoms, float res){ - - projectToTensor<<<1, num_atom_types>>>( coords, num_atoms_of_type, offsets, - volume, spatial_dim, res); - +/* +Input: + gradOutput: gradient of selected features + coords: coordinates in a flat array + num_atoms: number of atoms + spatial_dim: volume 3d array real size + res: volume 3d array resolution + +Output: + gradInput: 3d array from which we selected +*/ + int feature_index = threadIdx.x; + float *feature_volume = gradInput + feature_index * spatial_dim*spatial_dim*spatial_dim; + float *feature_grad = gradOutput + feature_index * max_num_atoms; + for(int atom_idx = 0; atom_idx=0)&&(y=0)&&(z=0)){ + uint idx = z + y*spatial_dim + x*spatial_dim*spatial_dim; + feature_volume[idx] = feature_grad[atom_idx]; + } + } } -template -void gpu_computeVolume2Coords( T *coords, - T* grad, - int *num_atoms_of_type, - int *offsets, - T *volume, - int spatial_dim, - int num_atom_types, - float res){ - projectFromTensor<<<1, num_atom_types>>>(coords, grad, num_atoms_of_type, offsets, - volume, spatial_dim, res); +void gpu_coordSelect( float *features, int num_features, + float* volume, int spatial_dim, + float *coords, int num_atoms, int max_num_atoms, + float res){ + + coordSelect<<<1, num_features>>>( features, + volume, spatial_dim, + coords, num_atoms, max_num_atoms, + res); } -void gpu_selectFromTensor( float *features, int num_features, - float* volume, int spatial_dim, +void gpu_coordSelectGrad( float *gradOutput, int num_features, + float* gradInput, int spatial_dim, float *coords, int num_atoms, int max_num_atoms, float res){ - selectFromTensor<<<1, num_features>>>( features, - volume, spatial_dim, + coordSelectGrad<<<1, num_features>>>( gradOutput, + gradInput, spatial_dim, coords, num_atoms, max_num_atoms, res); } +/* +template void gpu_computeVolume2Coords( float*, float*, int, float*, int, float); +template void gpu_computeVolume2Coords( double*, double*, int, double*, int, float); -template void gpu_computeVolume2Coords( float*, float*, int*, int*, float*, int, int, float); -template void gpu_computeVolume2Coords( double*, double*, int*, int*, double*, int, int, float); - -template void gpu_computeCoords2Volume(float*, int*, int*, float*, int, int, float); -template void gpu_computeCoords2Volume(double*, int*, int*, double*, int, int, float); \ No newline at end of file +template void gpu_computeCoords2Volume(float*, int, float*, int, float); +template void gpu_computeCoords2Volume(double*, int, double*, int, float); +*/ \ No newline at end of file diff --git a/Layers/Volume/Kernels.h b/Layers/Volume/Kernels.h index 0e1c537..d6c1ce7 100644 --- a/Layers/Volume/Kernels.h +++ b/Layers/Volume/Kernels.h @@ -1,25 +1,25 @@ -#include - +/* template void gpu_computeCoords2Volume( T *coords, - int *num_atoms_of_type, - int *offsets, + int num_atoms, T *volume, int spatial_dim, - int num_atom_types, float res); template void gpu_computeVolume2Coords( T *coords, T* grad, - int *num_atoms_of_type, - int *offsets, + int num_atoms, T *volume, int spatial_dim, - int num_atom_types, float res); +*/ +void gpu_coordSelect( float *features, int num_features, + float* volume, int spatial_dim, + float *coords, int num_atoms, int max_num_atoms, + float res); -void gpu_selectFromTensor( float *features, int num_features, - float* volume, int spatial_dim, +void gpu_coordSelectGrad( float *gradOutput, int num_features, + float* gradInput, int spatial_dim, float *coords, int num_atoms, int max_num_atoms, float res); \ No newline at end of file diff --git a/Layers/Volume/RotateGrid.cu b/Layers/Volume/RotateGrid.cu deleted file mode 100644 index a327e1a..0000000 --- a/Layers/Volume/RotateGrid.cu +++ /dev/null @@ -1,37 +0,0 @@ -#include "RotateGrid.h" - -__global__ void gpuRotatePoint(REAL *d_rotations, REAL *d_grid, int batch_size, int volume_size){ - int batch_id = blockIdx.x; - float *d_m = d_rotations + batch_id * 9; - - int volume = (volume_size*volume_size*volume_size); - int volume_slice = (volume_size*volume_size); - - int x = blockIdx.y; - int y = blockIdx.z; - int z = threadIdx.x; - - float d_v[3] = { 2.0*float(x + 0.5 - volume_size/2.0)/ float(volume_size), - 2.0*float(y + 0.5 - volume_size/2.0)/ float(volume_size), - 2.0*float(z + 0.5 - volume_size/2.0)/ float(volume_size) }; - - float *dst = d_grid + batch_id*volume*3 + x*volume_slice*3 + y*volume_size*3 + z*3; - - // x, y, z - // dst[0] = d_m[0]*d_v[0]+d_m[1]*d_v[1]+d_m[2]*d_v[2]; - // dst[1] = d_m[3]*d_v[0]+d_m[4]*d_v[1]+d_m[5]*d_v[2]; - // dst[2] = d_m[6]*d_v[0]+d_m[7]*d_v[1]+d_m[8]*d_v[2]; - - // z, y, x - dst[2] = d_m[0]*d_v[0]+d_m[1]*d_v[1]+d_m[2]*d_v[2]; - dst[1] = d_m[3]*d_v[0]+d_m[4]*d_v[1]+d_m[5]*d_v[2]; - dst[0] = d_m[6]*d_v[0]+d_m[7]*d_v[1]+d_m[8]*d_v[2]; - -} - -void cpu_RotateGrid(REAL *d_rotations, REAL *d_grid, int batch_size, int volume_size){ - - dim3 dim_special(batch_size, volume_size, volume_size); - gpuRotatePoint<<>>(d_rotations, d_grid, batch_size, volume_size); - -} \ No newline at end of file diff --git a/Layers/Volume/RotateGrid.h b/Layers/Volume/RotateGrid.h deleted file mode 100644 index 9334e2b..0000000 --- a/Layers/Volume/RotateGrid.h +++ /dev/null @@ -1,4 +0,0 @@ -#define REAL float -// #define REAL double - -void cpu_RotateGrid(REAL *d_rotations, REAL *d_grid, int batch_size, int volume_size); \ No newline at end of file diff --git a/Layers/Volume/Select/select_interface.cpp b/Layers/Volume/Select/select_interface.cpp index a07342c..ce9eefb 100644 --- a/Layers/Volume/Select/select_interface.cpp +++ b/Layers/Volume/Select/select_interface.cpp @@ -31,9 +31,43 @@ void SelectVolume_forward( torch::Tensor volume, torch::Tensor single_features = features[i]; int single_num_atoms = num_atoms[i].item().toInt(); - gpu_selectFromTensor( single_features.data(), num_features, - single_volume.data(), spatial_dim, - single_coords.data(), single_num_atoms, max_num_atoms, res); + gpu_coordSelect(single_features.data(), num_features, + single_volume.data(), spatial_dim, + single_coords.data(), single_num_atoms, max_num_atoms, res); } -} +} + + +void SelectVolume_backward( torch::Tensor gradOutput, + torch::Tensor gradInput, + torch::Tensor coords, + torch::Tensor num_atoms, + float res + ){ + CHECK_GPU_INPUT(gradOutput); + CHECK_GPU_INPUT(gradInput); + CHECK_GPU_INPUT(coords); + CHECK_GPU_INPUT_TYPE(num_atoms, torch::kInt); + if(coords.ndimension() != 2){ + ERROR("Incorrect input ndim"); + } + + + int batch_size = coords.size(0); + int num_features = gradOutput.size(1); + int spatial_dim = gradInput.size(2); + int max_num_atoms = coords.size(1)/3; + + #pragma omp parallel for + for(int i=0; i(), num_features, + single_gradInput.data(), spatial_dim, + single_coords.data(), single_num_atoms, max_num_atoms, res); + } +} diff --git a/Layers/Volume/Select/select_interface.h b/Layers/Volume/Select/select_interface.h index d6c0b15..dca7496 100644 --- a/Layers/Volume/Select/select_interface.h +++ b/Layers/Volume/Select/select_interface.h @@ -5,4 +5,11 @@ void SelectVolume_forward( torch::Tensor volume, torch::Tensor num_atoms, torch::Tensor features, float res + ); + +void SelectVolume_backward( torch::Tensor gradOutput, + torch::Tensor gradInput, + torch::Tensor coords, + torch::Tensor num_atoms, + float res ); \ No newline at end of file diff --git a/Layers/Volume/TypedCoords2Volume/typedcoords2volume_interface.cpp b/Layers/Volume/TypedCoords2Volume/typedcoords2volume_interface.cpp index 5e234a4..cc2d4e0 100644 --- a/Layers/Volume/TypedCoords2Volume/typedcoords2volume_interface.cpp +++ b/Layers/Volume/TypedCoords2Volume/typedcoords2volume_interface.cpp @@ -1,59 +1,83 @@ #include "typedcoords2volume_interface.h" #include #include -#include +#include #include +//#define _OPENMP +//#include + void TypedCoords2Volume_forward( torch::Tensor input_coords, torch::Tensor volume, - torch::Tensor num_atoms_of_type, - torch::Tensor offsets, - float resolution){ - int num_atom_types=11; + torch::Tensor num_atoms, + float resolution, int num_neighbours, + torch::Tensor gridParticleHash, + torch::Tensor gridParticleIndex, + torch::Tensor cellStart, + torch::Tensor cellStop, + torch::Tensor sortedPos){ CHECK_GPU_INPUT(volume); CHECK_GPU_INPUT(input_coords); - CHECK_GPU_INPUT_TYPE(num_atoms_of_type, torch::kInt); - CHECK_GPU_INPUT_TYPE(offsets, torch::kInt); + CHECK_CPU_INPUT_TYPE(num_atoms, torch::kInt); + + CHECK_GPU_INPUT_TYPE(gridParticleHash, torch::kInt64); + CHECK_GPU_INPUT_TYPE(gridParticleIndex, torch::kInt64); + CHECK_GPU_INPUT_TYPE(cellStart, torch::kInt64); + CHECK_GPU_INPUT_TYPE(cellStop, torch::kInt64); + CHECK_GPU_INPUT(sortedPos); + if(input_coords.ndimension() != 2){ ERROR("Incorrect input ndim"); } int batch_size = input_coords.size(0); - - #pragma omp parallel for - for(int i=0; i(); + // #pragma omp parallel for + // for(int i=0; i( single_input_coords.data(), - single_num_atoms_of_type.data(), - single_offsets.data(), - single_volume.data(), single_volume.size(1), num_atom_types, resolution); + a_num_atoms[i], + single_volume.data(), + single_volume.size(1), resolution, num_neighbours, + (single_particleHash[0]).data(), + (single_particleHash[1]).data(), + (single_particleIndex[0]).data(), + (single_particleIndex[1]).data(), + single_cellStart.data(), + single_cellStop.data(), + single_sortedPos.data(), + device); })); - } + }}); } void TypedCoords2Volume_backward( torch::Tensor grad_volume, torch::Tensor grad_coords, torch::Tensor coords, - torch::Tensor num_atoms_of_type, - torch::Tensor offsets, - float resolution){ - int num_atom_types=11; + torch::Tensor num_atoms, + float resolution, + int num_neighbours){ CHECK_GPU_INPUT(grad_volume); CHECK_GPU_INPUT(grad_coords); CHECK_GPU_INPUT(coords); - CHECK_GPU_INPUT_TYPE(num_atoms_of_type, torch::kInt); - CHECK_GPU_INPUT_TYPE(offsets, torch::kInt); + CHECK_CPU_INPUT_TYPE(num_atoms, torch::kInt); if(grad_coords.ndimension() != 2){ ERROR("Incorrect input ndim"); } int batch_size = grad_coords.size(0); + auto a_num_atoms = num_atoms.accessor(); #pragma omp parallel for for(int i=0; i(single_coords.data(), single_grad_coords.data(), - single_num_atoms_of_type.data(), - single_offsets.data(), + a_num_atoms[i], single_grad_volume.data(), - single_grad_volume.size(1), num_atom_types, resolution); + single_grad_volume.size(1), resolution, num_neighbours); })); } } + + diff --git a/Layers/Volume/TypedCoords2Volume/typedcoords2volume_interface.h b/Layers/Volume/TypedCoords2Volume/typedcoords2volume_interface.h index 199f14b..5facd1a 100644 --- a/Layers/Volume/TypedCoords2Volume/typedcoords2volume_interface.h +++ b/Layers/Volume/TypedCoords2Volume/typedcoords2volume_interface.h @@ -1,12 +1,17 @@ #include void TypedCoords2Volume_forward( torch::Tensor input_coords, - torch::Tensor volume, - torch::Tensor num_atoms_of_type, - torch::Tensor offsets, - float resolution); + torch::Tensor volume, + torch::Tensor num_atoms, + float resolution, int num_neighbours, + torch::Tensor gridParticleHash, + torch::Tensor gridParticleIndex, + torch::Tensor cellStart, + torch::Tensor cellStop, + torch::Tensor sortedPos); + void TypedCoords2Volume_backward( torch::Tensor grad_volume, - torch::Tensor grad_coords, - torch::Tensor coords, - torch::Tensor num_atoms_of_type, - torch::Tensor offsets, - float resolution); \ No newline at end of file + torch::Tensor grad_coords, + torch::Tensor coords, + torch::Tensor num_atoms, + float resolution, + int num_neighbours); \ No newline at end of file diff --git a/Layers/Volume/VolumeConv.cu b/Layers/Volume/VolumeConv.cu deleted file mode 100644 index 550cf28..0000000 --- a/Layers/Volume/VolumeConv.cu +++ /dev/null @@ -1,137 +0,0 @@ -#include "VolumeConv.h" -#include -#include - -#define gpuFFTErrchk(err) { __cufftSafeCall(err, __FILE__, __LINE__); } -#define gpuErrchk(err) { gpuAssert(err, __FILE__, __LINE__); } -inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=true) -{ - if (code != cudaSuccess) - { - fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line); - if (abort) exit(code); - } -} - -static const char *_cudaGetErrorEnum(cufftResult error) -{ - switch (error) - { - case CUFFT_SUCCESS: - return "CUFFT_SUCCESS"; - - case CUFFT_INVALID_PLAN: - return "CUFFT_INVALID_PLAN"; - - case CUFFT_ALLOC_FAILED: - return "CUFFT_ALLOC_FAILED"; - - case CUFFT_INVALID_TYPE: - return "CUFFT_INVALID_TYPE"; - - case CUFFT_INVALID_VALUE: - return "CUFFT_INVALID_VALUE"; - - case CUFFT_INTERNAL_ERROR: - return "CUFFT_INTERNAL_ERROR"; - - case CUFFT_EXEC_FAILED: - return "CUFFT_EXEC_FAILED"; - - case CUFFT_SETUP_FAILED: - return "CUFFT_SETUP_FAILED"; - - case CUFFT_INVALID_SIZE: - return "CUFFT_INVALID_SIZE"; - - case CUFFT_UNALIGNED_DATA: - return "CUFFT_UNALIGNED_DATA"; - } - - return ""; -} - -inline void __cufftSafeCall(cufftResult err, const char *file, const int line) -{ - if( CUFFT_SUCCESS != err) { - fprintf(stderr, "CUFFT error in file '%s', line %d\n %s\nerror %d: %s\nterminating!\n",__FILE__, __LINE__,err, \ - _cudaGetErrorEnum(err)); \ - cudaDeviceReset(); exit(1); \ - } -} - -#define WARP_SIZE 32 - -__global__ void conjMul(cufftComplex *c_volume1, cufftComplex *c_volume2, cufftComplex *c_output, int batch_size, int volume_size, bool conj){ - uint batch_idx = blockIdx.x; - uint warp_idx = blockIdx.y; - uint thread_idx = threadIdx.x; - int reduced_volume_size = (volume_size/2 + 1); - uint vol = volume_size*volume_size*reduced_volume_size; - uint full_vol = volume_size*volume_size*volume_size; - uint memory_idx = batch_idx*vol + warp_idx*WARP_SIZE + thread_idx; - - if( (warp_idx*WARP_SIZE + thread_idx) >= vol) //out of volume - return; - REAL re, im; - if(conj){ - re = c_volume1[memory_idx].x * c_volume2[memory_idx].x + c_volume1[memory_idx].y * c_volume2[memory_idx].y; - im = -c_volume1[memory_idx].x * c_volume2[memory_idx].y + c_volume1[memory_idx].y * c_volume2[memory_idx].x; - }else{ - re = c_volume1[memory_idx].x * c_volume2[memory_idx].x - c_volume1[memory_idx].y * c_volume2[memory_idx].y; - im = c_volume1[memory_idx].x * c_volume2[memory_idx].y + c_volume1[memory_idx].y * c_volume2[memory_idx].x; - } - // printf("%f %f \n", c_volume1[memory_idx].x, c_volume2[memory_idx].y); - c_output[memory_idx].x = re*2.0/full_vol; - c_output[memory_idx].y = im*2.0/full_vol; - -} - -void cpu_VolumeConv(REAL *d_volume1, REAL *d_volume2, REAL *d_output, int batch_size, int volume_size, bool conjugate){ - cufftHandle plan_fwd, plan_bwd; - cufftComplex *c_volume1, *c_volume2, *c_output; - int reduced_volume_size = (volume_size/2 + 1); - - gpuErrchk(cudaMalloc((void**)&c_volume1, sizeof(cufftComplex)*batch_size*volume_size*volume_size*reduced_volume_size)); - gpuErrchk(cudaMalloc((void**)&c_volume2, sizeof(cufftComplex)*batch_size*volume_size*volume_size*reduced_volume_size)); - gpuErrchk(cudaMalloc((void**)&c_output, sizeof(cufftComplex)*batch_size*volume_size*volume_size*volume_size)); - // printf("Finished malloc\n"); - int dimensions_real[] = {volume_size, volume_size, volume_size}; - int dimensions_complex[] = {volume_size, volume_size, reduced_volume_size}; - int batch_volume_real = volume_size*volume_size*volume_size; - int batch_volume_complex = volume_size*volume_size*reduced_volume_size; - int inembed[] = {volume_size, volume_size, volume_size}; - int onembed[] = {volume_size, volume_size, reduced_volume_size}; - // printf("Started plan\n"); - gpuFFTErrchk(cufftPlanMany(&plan_fwd, 3, dimensions_real, - inembed, 1, batch_volume_real, - onembed, 1, batch_volume_complex, - CUFFT_R2C, batch_size)); - - gpuFFTErrchk(cufftPlanMany( &plan_bwd, 3, dimensions_real, - onembed, 1, batch_volume_complex, - inembed, 1, batch_volume_real, - CUFFT_C2R, batch_size)); - // printf("Started fft\n"); - gpuFFTErrchk(cufftExecR2C(plan_fwd, d_volume1, c_volume1)); - gpuFFTErrchk(cufftExecR2C(plan_fwd, d_volume2, c_volume2)); - // printf("Finished fft\n"); - - dim3 dim_special(batch_size, batch_volume_complex/WARP_SIZE + 1); - conjMul<<>>(c_volume1, c_volume2, c_output, batch_size, volume_size, conjugate); - gpuErrchk( cudaPeekAtLastError() ); - gpuErrchk( cudaDeviceSynchronize() ); - - // printf("Finished mul\n"); - gpuFFTErrchk(cufftExecC2R(plan_bwd, c_output, d_output)); - // printf("Finished invfft\n"); - - gpuFFTErrchk(cufftDestroy(plan_fwd)); - gpuFFTErrchk(cufftDestroy(plan_bwd)); - gpuErrchk(cudaFree(c_volume1)); - gpuErrchk(cudaFree(c_volume2)); - gpuErrchk(cudaFree(c_output)); - // printf("Finished free\n"); -} - - diff --git a/Layers/Volume/VolumeConv.h b/Layers/Volume/VolumeConv.h deleted file mode 100644 index 416b6d3..0000000 --- a/Layers/Volume/VolumeConv.h +++ /dev/null @@ -1,7 +0,0 @@ - -#define REAL float -// #define REAL double - - -void cpu_VolumeConv(REAL *d_volume1, REAL *d_volume2, REAL *d_output, int batch_size, int volume_size, bool conjugate); -// Computes F(tau) = \int vol1(r)vol2(tau-r) dr diff --git a/Layers/Volume/VolumeConvolution/volumeConvolution_interface.cpp b/Layers/Volume/VolumeConvolution/volumeConvolution_interface.cpp deleted file mode 100644 index f153097..0000000 --- a/Layers/Volume/VolumeConvolution/volumeConvolution_interface.cpp +++ /dev/null @@ -1,52 +0,0 @@ -#include "volumeConvolution_interface.h" -#include -#include -#include - - -void VolumeConvolution_forward( torch::Tensor volume1, - torch::Tensor volume2, - torch::Tensor output){ - CHECK_GPU_INPUT_TYPE(volume1, torch::kFloat); - CHECK_GPU_INPUT_TYPE(volume2, torch::kFloat); - CHECK_GPU_INPUT_TYPE(output, torch::kFloat); - if(volume1.ndimension()!=4){ - ERROR("incorrect input dimension"); - } - cpu_VolumeConv( volume1.data(), - volume2.data(), - output.data(), - volume1.size(0), - volume1.size(1), - true); -} -void VolumeConvolution_backward( torch::Tensor gradOutput, - torch::Tensor gradVolume1, - torch::Tensor gradVolume2, - torch::Tensor volume1, - torch::Tensor volume2){ - CHECK_GPU_INPUT_TYPE(gradVolume1, torch::kFloat); - CHECK_GPU_INPUT_TYPE(gradVolume2, torch::kFloat); - CHECK_GPU_INPUT_TYPE(gradOutput, torch::kFloat); - CHECK_GPU_INPUT_TYPE(volume1, torch::kFloat); - CHECK_GPU_INPUT_TYPE(volume2, torch::kFloat); - if(gradOutput.ndimension()!=4){ - ERROR("incorrect input dimension"); - } - - cpu_VolumeConv( gradOutput.data(), - volume2.data(), - gradVolume1.data(), - volume1.size(0), - volume1.size(1), - false); - - cpu_VolumeConv( volume1.data(), - gradOutput.data(), - gradVolume2.data(), - volume1.size(0), - volume1.size(1), - true); -} - - diff --git a/Layers/Volume/VolumeConvolution/volumeConvolution_interface.h b/Layers/Volume/VolumeConvolution/volumeConvolution_interface.h deleted file mode 100644 index 4c2404b..0000000 --- a/Layers/Volume/VolumeConvolution/volumeConvolution_interface.h +++ /dev/null @@ -1,9 +0,0 @@ -#include -void VolumeConvolution_forward( torch::Tensor volume1, - torch::Tensor volume2, - torch::Tensor output); -void VolumeConvolution_backward( torch::Tensor gradOutput, - torch::Tensor gradVolume1, - torch::Tensor gradVolume2, - torch::Tensor volume1, - torch::Tensor volume2); \ No newline at end of file diff --git a/Layers/Volume/VolumeRotation/volumeRotation_interface.cpp b/Layers/Volume/VolumeRotation/volumeRotation_interface.cpp deleted file mode 100644 index d351d1f..0000000 --- a/Layers/Volume/VolumeRotation/volumeRotation_interface.cpp +++ /dev/null @@ -1,17 +0,0 @@ -#include "volumeRotation_interface.h" -#include -#include -#include - - -void VolumeGenGrid( torch::Tensor rotations, torch::Tensor grid){ - CHECK_GPU_INPUT_TYPE(rotations, torch::kFloat32); - CHECK_GPU_INPUT_TYPE(grid, torch::kFloat32); - if(rotations.ndimension()!=3 || grid.ndimension()!=5){ - ERROR("incorrect input dimension"); - } - int batch_size = rotations.size(0); - int size = grid.size(1); - cpu_RotateGrid(rotations.data(), grid.data(), batch_size, size); -} - diff --git a/Layers/Volume/VolumeRotation/volumeRotation_interface.h b/Layers/Volume/VolumeRotation/volumeRotation_interface.h deleted file mode 100644 index 3b91d03..0000000 --- a/Layers/Volume/VolumeRotation/volumeRotation_interface.h +++ /dev/null @@ -1,2 +0,0 @@ -#include -void VolumeGenGrid( torch::Tensor rotations, torch::Tensor grid); \ No newline at end of file diff --git a/Layers/Volume/main.cpp b/Layers/Volume/main.cpp index 01330a6..1917c77 100644 --- a/Layers/Volume/main.cpp +++ b/Layers/Volume/main.cpp @@ -2,8 +2,6 @@ #include #include #include