Format_pdb.cpp

#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <numeric>
#include <cmath>
#include <algorithm>
#include <list>
#include <map>
#include <iterator>

using namespace std;

int main(int argc, char **argv)
{
	if(argc!=4)
	{
		cout<<"Incorrect Useage:"<<std::endl;
		cout<<"Correct useage is: "<<argv[0]<<" input.pdb  original.pdb output.pdb"<<endl;
		return 1;
	}
	std::ifstream file(argv[2]);
	std::stringstream pdb;
	if(file.is_open())
	{
		pdb << file.rdbuf();
		file.close();
	}
	else
	{
		cout<<"Unable to open "<<argv[2]<<endl;
		return 1;		
	}
	std::string line;
	vector<string> atoms;
	while(std::getline(pdb,line))	if(line.substr(0,4) == "ATOM")	atoms.push_back(line);
	file.open(argv[1]);
	std::stringstream data;
	if(file.is_open())
	{
		data << file.rdbuf();
		file.close();
	}
	else
	{
		cout<<"Unable to open "<<argv[1]<<endl;
		return 1;
	}
	
	vector<string> new_atoms;
	vector<string> water;
	while(std::getline(data,line))
	{
		if(line.substr(0,4) == "ATOM")
		{
			if(line.substr(17,3) == "HIE") line.replace(17,3, "HIS"); 
			if(line.substr(17,3) == "HID") line.replace(17,3, "HIS");
			if(line.substr(17,3) == "HIP") line.replace(17,3, "HIS");
			new_atoms.push_back(line);
		}
		if(line.substr(0,6) == "HETATM" && line.substr(17,3) == "WAT")
		{
			line.replace(17,3, "HOH");
			water.push_back(line);
		}
	}
	
	std::vector <string> res_list1, res_list2;
	std::vector<string>::iterator it;
	for (size_t i=0; i<atoms.size();++i)	res_list1.push_back(atoms.at(i).substr(17,9));
	it = std::unique (res_list1.begin(), res_list1.end());
	res_list1.resize( std::distance(res_list1.begin(),it) );
	
	for (size_t i=0; i<new_atoms.size();++i)	res_list2.push_back(new_atoms.at(i).substr(17,9));
	it = std::unique (res_list2.begin(), res_list2.end());
	res_list2.resize( std::distance(res_list2.begin(),it) );
	
	map <string,string> map1;
		
	for (size_t i=0; i<res_list2.size();++i)	map1[res_list2[i].substr(5,4)] = res_list1[i].substr(5,4);
	for (size_t i=0; i<new_atoms.size();++i)	new_atoms[i]=new_atoms.at(i).replace(22,4,map1[new_atoms.at(i).substr(22,4)]);
	
	vector<string> water_needed;
	size_t atom_count= new_atoms.size()+2;
	for (size_t i=0; i<water.size();++i)
	{
		auto x1=stof(water.at(i).substr(30,8));
		auto y1=stof(water.at(i).substr(38,8));
		auto z1=stof(water.at(i).substr(46,8));
		bool is_needed = false;
		#pragma omp parallel for
		for (size_t j=0; j<new_atoms.size();++j)
		{
			auto x2=stof(new_atoms.at(j).substr(30,8));
			auto y2=stof(new_atoms.at(j).substr(38,8));
			auto z2=stof(new_atoms.at(j).substr(46,8));
			if(!is_needed)
			{
				auto d = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
				if (d <= 4.50)
				{
					is_needed = true;
				}
			}
		}
		if(is_needed == true)
		{	atom_count++;
			string k=to_string(atom_count);
			while(k.length()<5) k.insert(0,1,' ');
			water_needed.push_back(water.at(i).replace(6,5,k));
		}
	}
	
	std::ofstream outFile;
    outFile.open(argv[3]);
	for (const auto &e : new_atoms) outFile << e << "\n";
	for (const auto &e : water_needed) outFile << e << "\n";
	
	
	return 0;
}