#!/usr/bin/perl

# 23rd Oct 2015 Pratiti Bhadra
# INPUT: Protein ID (PDB file name) 
# OUTPUT :  .itp file and .top file as per CGMM force field
# $ perl pdb2ItpGroTop_multimer_multidomain.pl XXXX

use warnings; 
use strict;

my $filename = $ARGV[0]; # input only name PDB XXXX
my $pdbfile = $ARGV[0]. '.pdb';
my $trac_res = $filename.'_res_trac.txt';
my @chain;
my $len_word;
my @chainid;
my @ChEnd;
my @ChStart;
my @domainStart;
my @domainEnd;


##########################
ReadArgument();



##########################
sub ReadArgument{

    if ($ARGV[0] eq "") {
	die "\nPlease enter Protein PDB ID as argument \nSyntax: pdb2ItpGroTop_multimer_multidomain.pl proteinPDBID\n\n";
    	
    }
    elsif ($ARGV[0] eq "-help") {
    	&PrintHelp;
    }
    elsif($ARGV[0] eq "-h")
    {
	die"Syntax: pdb2ItpGroTop_multimer_multidomain.pl proteinPDBID\n";
    }
    else
    {

	my $filelen = length($filename);
	if($filelen != 4)
	{
		print"\n## Please enter a PDB ID of four character\n\n";
		exit;
	}	
	my $pdbfile = $ARGV[0]. '.pdb';

	my $filepath = './'.$pdbfile;
	unless (-e $filepath) {
 	print "\n## File Doesn't Exist!\n\n";
	exit;
 	} 
	
	
	
	
	my $ispdb = 0;	
	
	my @chainidX;
	my $len_wordX;
	my @wordX;

	open ( FXPDB, "< ./$pdbfile") or die "cannot open the file $!";
	while(<FXPDB>)
	{
		if(/^COMPND/)
		{ 
			@wordX = split(/\s+/,$_);
			
					
			if($wordX[2] =~ /CHAIN:/ )
			{
				$len_wordX = scalar(@wordX);
				
				for(my $i = 3; $i < $len_wordX; $i++)
				{
					my $tempX = substr($wordX[$i],0,1,);
							
					push @chainidX, $tempX;
				}
				
			}
			
			
		}
	}
	close FXPDB;
	for(my $chcnt= 0; $chcnt<scalar(@chainidX); $chcnt++) # for each chain check residue gap
	{
		open ( FXPDB, "< ./$pdbfile") or die "cannot open the file $!";
		my @res_array;
		while(<FXPDB>)
		{

			if (/^ATOM/)
         		{
				if(substr($_,13,2) eq 'CA' && substr($_,21,1) eq $chainidX[$chcnt])
				{
					$ispdb = 1;
					push @res_array,sprintf("%d",substr($_,22,4));
		
				}
			}
			if(/^ENDMDL/)
			{
				last;
			}
		}
		close FXPDB;
		
				
		if((scalar(@res_array)-1) != $res_array[scalar(@res_array)-1]-$res_array[0])
		{
			die"\nThere is Residue gap or Alternate location indicator in PDB file in chain $chainidX[$chcnt]\n\n";
		}


	}

		
	
	if($ispdb == 0)
	{
		die"\nThere is no information about Calpha ATOM\n\n";
	}
	


	chain_number();

	CGMM_res_no_exclude_bond_multichian();

	linker_region();

	itptopgro();
    }
    
}

sub PrintHelp {
  open HELP, "| more";
  print <<EndHelp;
NAME
        pdb2ItpGroTop_multimer_multidomain.pl

DESCRIPTION
	Create topology (.itp, .top) and coordinate (.gro) file of monomer single domain protein for simulation run 
	http://pallab.serc.iisc.ernet.in/CGMM/

AUTHOR
	Pratiti Bhadra (pratiti.bhadra[@]gmail.com)  

OPTIONS
		       
	-h 
		Give syntax. This argument must be the first.

	-help 
		Give detailed help (print this message). This argument 
		must be the first.

EXAMPLE
        perl pdb2ItpGroTop_multimer_multidomain.pl 1TRL
	perl pdb2ItpGroTop_multimer_multidomain.pl 1D1D

EndHelp

close HELP;
die "\n";

}  #  End of PrintHelp
sub chain_number {
	
	open ( FPDB, "< $pdbfile") or die "cannot open the file $!";
	
	open ( FTR, "> ./$trac_res") or die "cannot open the file $!";
	

	my $count = 0;
	
	while(<FPDB>)	
	{
		my @word;
		chomp $_; 
					
		if(/^COMPND/)
		{ 
			@word = split(/\s+/,$_);
			
					
			if($word[2] =~ /CHAIN:/ )
			{
				$len_word = scalar(@word);
				
				for(my $i = 3; $i < $len_word; $i++)
				{
					my $temp = substr($word[$i],0,1,);
							
					push @chainid, $temp;
				}
				
			}
			
			
		}
		
		if(/^ATOM/)
		{
			my @col;			
			if(substr($_,13,2) eq 'CA' )
                 	{
				@col = split(/\s+/,$_);
				$count++;

				printf FTR "%5d %s %s\n", $count,$col[4],$col[5];
				
			}
		}
		if(/^ENDMDL/)
		{
			last;
		}
	}
	close FPDB;
	close FTR;
			
}

sub CGMM_res_no_exclude_bond_multichian {

	open ( F2, "> ./exclude_res.txt") or die "cannot open the file $!";
	
	
	for(my $k=1; $k<scalar(@chainid); $k++)
	{
		my $flag = 0;
				
		open ( FTR, "< ./$trac_res") or die "cannot open the file $!";	
		while(<FTR>)
		{
			chomp $_;
			
			my $ch_comp = substr($_,6,1); 
			my $no = substr($_,0,5);			

			if($ch_comp eq $chainid[$k])
			{
				$flag = 1;
				my $a = sprintf("%d",$no);
				
				printf F2 "%5d%5d\n", ($a-1), $a;
				push @ChEnd, ($a-1); # end of chain
				push @ChStart, $a; # start of chain
			}
			if($flag == 1)
			{
				$flag = 0;				
				last;
				
			}

		}
		close FTR;
			 
	}

	close F2;
  

}

sub linker_region {
	
		
	print"\n\n#########\n\nEnter linker regions residue numbers (PDB file) (NNN NNN#NNN NNN format, for each chain); \nExample: \n\nIf chain A has one linker region \n(linker region residue number from 78 to 86)\nChain A:\n78 86. \n\nIf chain A has more that one linker region (ther is three linker regions residue number from 78 to 86, from 105 to 118 and from 130 to 140)\nChain A:\n78 86#105 118#130 140\n\nIf there is no linker region \nChain A:\n0 0\n\n#########\n\nPlease provide your input (domain linker region of each chain):\n";	
		
	for(my $cnt=0; $cnt<scalar(@chainid); $cnt++)
	{
		print"Chain $chainid[$cnt]:\n";
		my $age = <STDIN>; 
		my @temp_arr = split(/#/,$age);
                my @temp_arr2;
		
		for(my $c = 0; $c<scalar(@temp_arr); $c++)
		{ 
			@temp_arr2 = split(/\s+/,$temp_arr[$c]);
					
			if(sprintf("%d",$temp_arr2[0]) <= sprintf("%d",$temp_arr2[1]))
			{
				
				open(FTRAC, "<./$trac_res");
				while(my $trc = <FTRAC>)
				{
					chomp($trc);
					my @line = split(/\s+/,$trc);
					my $line3_no = sprintf("%d",$line[3]);
					
					if(($chainid[$cnt] eq $line[2])&& ($line3_no == $temp_arr2[0]) )
					{
						push @domainStart ,sprintf("%d", $line[1]);
					}
					if(($chainid[$cnt] eq $line[2])&& ($line3_no == $temp_arr2[1]) )
					{
						push @domainEnd ,sprintf("%d", $line[1]);
					}

				}
				close FTRAC;
				
			}
			else
			{
				print"Error! in linker region!!!\n";
				exit;
		
			}
		}

	
	}
}

sub itptopgro {

	my $a = scalar(@ChStart);

	# check is it metalloprotein
	system("grep -A50 'M RES CSSEQI METAL' $pdbfile | grep 'REMARK 620' > ./test.txt");
	open ( FXX, "< ./test.txt") or die "cannot open the file $!";
	
	my $metal_flg = 0;
	my $metal_atom;
	while(<FXX>)
	{
		chomp($_);
		my @array = split(/\s+/,$_);
		for(my $r=0;$r<scalar(@array); $r++)
		{
			if($array[$r] eq "CA")
			{
				$metal_flg = 1;
				$metal_atom = 'CA';
			}
			if($array[$r] eq "CU")
			{
				$metal_flg = 1;
				$metal_atom = 'CU';
			}
			if($array[$r] eq "ZN")
			{
				$metal_flg = 1;
				$metal_atom = 'ZN';
			}
			if($array[$r] eq "FE")
			{
				$metal_flg = 1;
				$metal_atom = 'FE';
			}
			if($array[$r] eq "MG")
			{
				$metal_flg = 1;
				$metal_atom = 'MG';
			}
			
		}

	}
	if($metal_flg == 0)
	{
		print"\nThere is no Metal ion. Please check.\n\n";

	}
	close FXX;
	system("rm ./test.txt");


	system "./secstr.exc ./$pdbfile"; 

	print("\n\nNAME
		pdb2ItpGroTop_multimer_multidomain.pl

	DESCRIPTION
		Create topology (.itp, .top) and coordinate (.gro) file of monomer single domain protein for simulation run 
		http://pallab.serc.iisc.ernet.in/CGMM/

	AUTHOR
		Pratiti Bhadra (pratiti.bhadra[@]gmail.com)  

	OPTIONS
			       
		-h 
			Give syntax. This argument must be the first.

		-help 
			Give detailed help (print this message). This argument 
			must be the first.

	EXAMPLE
		perl pdb2ItpGroTop_multimer_multidomain.pl 1TRL
		perl pdb2ItpGroTop_multimer_multidomain.pl 1D1D\n\n");
	print("\n!! WAIT !!\n");

	my $rinfile = $filename.'.rin';
	my $rinfile1 = $filename.'_001.rin';
	if (-e $rinfile1){ system("mv $rinfile1 $rinfile")};


	# secondary structure
	my $secstrfile = $filename . '.ssd';
	# residue one letter code
	my $res_1code_file = $filename . '.seq'; 
	# atom type file
	my $atom_type_file = $filename .'_atom_type.txt';;
	# residue three letter code
	my $res_3code_file = 'amino_acid_3letter_code.txt';

	open ( FH, "< $rinfile") or die "cannot open the file $!";
	open (F3, "> $secstrfile")or die "cannot open the file $!";
	open (F4, "> $atom_type_file")or die "cannot open the file $!";
	open (F5, "> $res_3code_file")or die "cannot open the file $!";
	open (F6, "> $res_1code_file")or die "cannot open the file $!";

	open (F, "> c.txt")or die "cannot open the file $!";

	my $list;
	my $ch_id;
	my $d;
	my @secstr;
	my @atom;
	my @bond_type;
	my @ori_resid_no;
	my @cgmm_resid_no;
	my @ori_chid;
	my $ss;

	while(<FH>)
	{
		# secondary structure single letter
		my $ss=uc(substr($_,14,1));
		my $rno1 = sprintf("%d",substr($_,0,4));
		my $rno2 = sprintf("%d",substr($_,9,4));
		my $rno3 = substr($_,8,1);
		push @ori_resid_no, $rno2;
		push @cgmm_resid_no, $rno1;
		push @ori_chid, $rno3;
		 
		if ($ss eq " ")
		{     $ss = "L";  }
		if($ss eq "B")
		{     $ss = "E" }
		$list = $list.$ss;
		print F3 "$ss"; 
		push @secstr, $ss;
		$ch_id = substr($_,8,1);
		# one letter amino acid code file and three letter amino acid code file
		my $resi=substr($_,4,3);
		print F5 "$resi\n";
		if($resi eq "ALA")
        	{   $d = "A"; }
        	if($resi eq "ARG")
        	{   $d = "R"; }
        	if($resi eq "ASN")
        	{   $d = "N"; }
        	if($resi eq "ASP")
        	{   $d = "D"; }
        	if($resi eq "CYS" || $resi eq "CSO" || $resi eq "CSX" || $resi eq "CME")
        	{   $d = "C"; }
        	if($resi eq "GLN")
        	{   $d = "Q"; }
        	if($resi eq "GLU" || $resi eq "PCA")
        	{   $d = "E"; }
        	if($resi eq "GLY")
        	{   $d = "G"; }
        	if($resi eq "HIS")
        	{   $d = "H"; }
        	if($resi eq "ILE")
        	{   $d = "I"; }
        	if($resi eq "LEU")
        	{   $d = "L"; }
        	if($resi eq "LYS" || $resi eq "MCL" || $resi eq "LLP" || $resi eq "KCX" || $resi eq "KPI")
        	{   $d = "K"; }
        	if($resi eq "MET" || $resi eq "MSE" || $resi eq "FME")
        	{  $d = "M"; }
        	if($resi eq "PHE")
        	{   $d = "F"; }
        	if($resi eq "PRO")
        	{   $d = "P"; }
        	if($resi eq "SER" || $resi eq "SEP" )
        	{   $d = "S"; }
        	if($resi eq "THR" || $resi eq "TPO" )
        	{   $d = "T"; }
        	if($resi eq "TRP")
        	{   $d = "W"; }
        	if($resi eq "TYR")
        	{   $d = "Y"; }
        	if($resi eq "VAL")
        	{   $d = "V"; }
		print F6 "$d\n";

		# atom type 
		my $tmp = 'C'.$d.$ss;
		print F4 "$tmp\n";
		 
		push @atom, $d;

	       # bond type
	       my $bd_tp = abs(sprintf("%8.3f",substr($_,29,7)));
	       printf F "$bd_tp\n";
	       if($bd_tp < 5)
	       {  push @bond_type, 0;  }
	       else
	       {  push @bond_type, 1; }         
		
	      

	}
	close F;
	# number of residue
	my $totalnores = length($list);

	close FH;
	close F3;
	close F4;
	close F5;
	close F6;

	# create triplet of protein
	open (F1, "> triplet.txt")or die "cannot open the file $!";
	my $i = 0;
	my $lx = scalar (@secstr);
	for( $i = 0;$i<$lx-2 ;$i++)
	{
	      	my $tempsrt = $secstr[$i].$secstr[$i+1].$secstr[$i+2];
		{
			my $i1 = sprintf("%4d",$i+1);
		        my $i2 = sprintf("%4d",$i+2); 
		        my $i3 = sprintf("%4d",$i+3);
			printf F1 "$i1 $i2 $i3 $tempsrt\n"  
		}
		           
	}
	close F1;
	# create quartet of protein
	open (F2, "> quartet.txt")or die "cannot open the file $!";
	$i = 0;
	for( $i = 0;$i<$lx-3 ;$i++)
	{
	      	my $tempsrt = $secstr[$i].$secstr[$i+1].$secstr[$i+2].$secstr[$i+3];
		{
			 my $i1 = sprintf("%4d",$i+1);
		         my $i2 = sprintf("%4d",$i+2); 
		         my $i3 = sprintf("%4d",$i+3);
		         my $i4 = sprintf("%4d",$i+4);
		         printf F2 "$i1 $i2 $i3 $i4 $tempsrt\n"  
		}
		           
	}
	close F2;	
	# create nonbonded pair
	open (F3, "> nonbonded_pair.txt")or die "cannot open the file $!";

	my $j;
	my $atom_number = scalar(@atom);

	for( $i = 1;$i<($atom_number-3);$i++)
	{
		    	
		for( $j = $i+2;$j<$atom_number-1;$j++)
	      	{
			my $tempsrt = $atom[$i].$atom[$j];
		        my $i1 = sprintf("%5d",$i+1);
		        my $j1 = sprintf("%5d",$j+1);
			printf F3 "$i1$j1 $tempsrt\n";              
		}
	
	}
	close F3;

	# metallo information 
	
	
	my @bsr_resno;
	my @bsr_resnm;
	my @bsr_atnm;
	my @bsr_chid;
	my @ori_at;
	

	my @cu_resno;
	my @Cux;
	my @Cuy;
	my @Cuz;
	my @Cun;
	my @Cuchid;
	my @cu_chid;

        if($metal_flg == 1 && $metal_atom eq 'CU')
        {
            my $fflag = 0;

            open ( FPDB, " < ./$pdbfile") or die "cannot open the file $!";
            while(my $fpdb = <FPDB>)
            {
                    chomp $fpdb;
                    if($fpdb =~ /^LINK/)
                    {
                            my @array_link = split(/\s+/,$fpdb);
			    if(($array_link[2] eq 'CU1' || $array_link[2] eq 'CU')  && $array_link[6] ne 'HOH')
                            {
                                    push @cu_resno,$array_link[4];
				    push @cu_chid,$array_link[3];
                                    push @bsr_resno,$array_link[8];
                                    push @bsr_resnm,$array_link[6];
				    push @bsr_chid,$array_link[7];
				    
                                    
                                    $fflag = 1;

                            }
                            if($array_link[2] ne 'HOH' && ($array_link[6] eq 'CU' || $array_link[6] eq 'CU1')  )
                            {
                                    push @cu_resno,$array_link[8];
				    push @cu_chid,$array_link[7];
                                    push @bsr_resno,$array_link[4];
                                    push @bsr_resnm,$array_link[2];
				    push @bsr_chid,$array_link[3];
                                    $fflag = 1;

                            }
                            
                    }
                    
                    if($fpdb =~ /^HETATM/)
                    {
                            
                            if( (substr($fpdb,17,3) eq ' CU' || substr($fpdb,17,3) eq 'CU1'))
                            {
                                    push @Cux,0.1*sprintf("%f",substr($fpdb,30,8));
                                    push @Cuy,0.1*sprintf("%f",substr($fpdb,38,8));
                                    push @Cuz,0.1*sprintf("%f",substr($fpdb,46,8));
                                    push @Cun,sprintf("%d",substr($fpdb,22,4));
				    push @Cuchid,substr($fpdb,21,1);
                                    
                                            
                            }
                                    
                    }
                    if($fpdb =~ /^ENDMDL/)
                    {
                                    last;
                    }

                    
            }
            close FPDB;
            
        }

	my @zn_resno;
	my @Znx;
	my @Zny;
	my @Znz;
	my @Znn;
	my @Znchid;
	my @zn_chid;
	if($metal_flg == 1 && $metal_atom eq 'ZN')
	{
		my $fflag = 0;

		
		open ( FPDB, " < $pdbfile") or die "cannot open the file $!";
		while(my $fpdb = <FPDB>)
		{
			chomp $fpdb;
			if($fpdb =~ /^LINK/)
			{
				my @array_link = split(/\s+/,$fpdb);
				if($array_link[2] eq 'ZN' && $array_link[6] ne 'HOH' )
				{
					push @zn_resno,$array_link[4];
					push @zn_chid,$array_link[3];
                                        push @bsr_resno,$array_link[8];
                                        push @bsr_resnm,$array_link[6];
				        push @bsr_chid,$array_link[7];
					
					
					$fflag = 1;

				}
				
				if($array_link[2] ne 'HOH' && $array_link[6] eq 'ZN' )
				{
					push @zn_resno,$array_link[8];
					push @zn_chid,$array_link[7];
                                        push @bsr_resno,$array_link[4];
                                        push @bsr_resnm,$array_link[2];
				        push @bsr_chid,$array_link[3];
					
					$fflag = 1;

				}
		
			}
	
			if($fpdb =~ /^HETATM/)
			{
		
				if(substr($fpdb,17,3) eq ' ZN')
				{
					push @Znx,0.1*sprintf("%f",substr($fpdb,30,8));
					push @Zny,0.1*sprintf("%f",substr($fpdb,38,8));
					push @Znz,0.1*sprintf("%f",substr($fpdb,46,8));
					push @Znn,sprintf("%d",substr($fpdb,22,4));
					push @Znchid,substr($fpdb,21,1);
			
				
				}
			
			}
			if($fpdb =~ /^ENDMDL/)
			{
					last;
			}

	
		}
		close FPDB;
	}

	my @fe_resno;
	my @Fex;
	my @Fey;
	my @Fez;
	my @Fen;
	my @Fechid;
	my @fe_chid;

	if($metal_flg == 1 && $metal_atom eq 'FE')
	{
		my $fflag = 0;

		open ( FPDB, " < $pdbfile") or die "cannot open the file $!";
		while(my $fpdb = <FPDB>)
		{
			chomp $fpdb;
			if($fpdb =~ /^LINK/)
			{
				my @array_link = split(/\s+/,$fpdb);
				if(($array_link[2] eq 'FE2' || $array_link[2] eq 'FE') && $array_link[6] ne 'HOH' )
				{
					push @fe_resno,$array_link[4];
					push @fe_chid,$array_link[3];
					push @bsr_resno,$array_link[8];
                                        push @bsr_resnm,$array_link[6];
				        push @bsr_chid,$array_link[7];
			
					$fflag = 1;

				}
				if($array_link[2] ne 'HOH' && ($array_link[6] eq 'FE2' || $array_link[6] eq 'FE'))
				{
					push @fe_resno,$array_link[8];
					push @fe_chid,$array_link[7];
					push @bsr_resno,$array_link[4];
                                        push @bsr_resnm,$array_link[2];
				        push @bsr_chid,$array_link[3];
					$fflag = 1;

				}
		
			}
	
			if($fpdb =~ /^HETATM/)
			{
		
				if( (substr($fpdb,17,3) eq 'FE2'|| substr($fpdb,17,3) eq ' FE'))
				{
					push @Fex,0.1*sprintf("%f",substr($fpdb,30,8));
					push @Fey,0.1*sprintf("%f",substr($fpdb,38,8));
					push @Fez,0.1*sprintf("%f",substr($fpdb,46,8));
					push @Fen,sprintf("%d",substr($fpdb,22,4));
					push @Fechid,substr($fpdb,21,1);
			
				
				}
			
			}
			if($fpdb =~ /^ENDMDL/)
			{
					last;
			}

	
		}
		close FPDB;
	}

	my @mg_resno;
	my @Mgx;
	my @Mgy;
	my @Mgz;
	my @Mgn;
	my @Mgchid;
	my @mg_chid;

	if($metal_flg == 1 && $metal_atom eq 'MG')
	{

		my $fflag = 0;

		open ( FPDB, " < $pdbfile") or die "cannot open the file $!";
		while(my $fpdb = <FPDB>)
		{
			chomp $fpdb;
			if($fpdb =~ /^LINK/)
			{
				my @array_link = split(/\s+/,$fpdb);
				if($array_link[2] eq 'MG' && $array_link[6] ne 'HOH')
				{
					push @mg_resno,$array_link[4];
					push @mg_chid,$array_link[3];
                                        push @bsr_resno,$array_link[8];
                                        push @bsr_resnm,$array_link[6];
				        push @bsr_chid,$array_link[7];
					my $sch_at = $array_link[5];
					push @ori_at,$sch_at;
					if($sch_at eq 'O' || $sch_at eq 'C' || $sch_at eq 'N' || $sch_at eq 'CA' )
					{
						push @bsr_atnm,'MC';
				
					}
					else
					{
						push @bsr_atnm,'SC';
				
					}
			
					$fflag = 1;

				}
				if($array_link[2] ne 'HOH' && $array_link[6] eq 'MG' )
				{
					push @mg_resno,$array_link[8];
					push @mg_chid,$array_link[7];
                                        push @bsr_resno,$array_link[4];
                                        push @bsr_resnm,$array_link[2];
				        push @bsr_chid,$array_link[3];
				        my $sch_at = $array_link[1];
					push @ori_at,$sch_at;
					if($sch_at eq 'O' || $sch_at eq 'C' || $sch_at eq 'N' || $sch_at eq 'MG' )
					{
						push @bsr_atnm,'MC';
					}
					else
					{
						push @bsr_atnm,'SC';
					}
					$fflag = 1;
			
				}
		
			}
	
			if($fpdb =~ /^HETATM/)
			{
		
				if( substr($fpdb,17,3) eq ' MG')
				{
					push @Mgx,0.1*sprintf("%f",substr($fpdb,30,8));
					push @Mgy,0.1*sprintf("%f",substr($fpdb,38,8));
					push @Mgz,0.1*sprintf("%f",substr($fpdb,46,8));
					push @Mgn,sprintf("%d",substr($fpdb,22,4));
					push @Mgchid,substr($fpdb,21,1);
			
				
				}
			
			}
			if($fpdb =~ /^ENDMDL/)
			{
					last;
			}

	
		}
		close FPDB;
	}

	my @ca_resno;
	my @Cax;
	my @Cay;
	my @Caz;
	my @Can;
	my @Cachid;
	my @ca_chid;
	
	if ($metal_flg == 1 && $metal_atom eq 'CA')
	{
		my $fflag = 0;

		open ( FPDB, " < $pdbfile") or die "cannot open the file $!";
		while(my $fpdb = <FPDB>)
		{
			chomp $fpdb;
	
			if($fpdb =~ /^LINK/)
			{
			    my @array_link = split(/\s+/,$fpdb);
			    
		     	    if($array_link[2] eq 'CA' && $array_link[6] ne 'HOH')
				{
					push @ca_resno,$array_link[4];
					push @ca_chid,$array_link[3];
                                        push @bsr_resno,$array_link[8];
                                        push @bsr_resnm,$array_link[6];
				        push @bsr_chid,$array_link[7];
					my $sch_at = $array_link[5];
					push @ori_at,$sch_at;
					if($sch_at eq 'O' || $sch_at eq 'C' || $sch_at eq 'N' || $sch_at eq 'CA' )
					{
						push @bsr_atnm,'MC';
				
					}
					else
					{
						push @bsr_atnm,'SC';
				
					}
			
					$fflag = 1;

				}
		            if($array_link[2] ne 'HOH' && $array_link[6] eq 'CA' )
				{
					push @ca_resno,$array_link[8];
					push @ca_chid,$array_link[7];
                                        push @bsr_resno,$array_link[4];
                                        push @bsr_resnm,$array_link[2];
				        push @bsr_chid,$array_link[3];
				        my $sch_at = $array_link[1];
					push @ori_at,$sch_at;
					if($sch_at eq 'O' || $sch_at eq 'C' || $sch_at eq 'N' || $sch_at eq 'CA' )
					{
						push @bsr_atnm,'MC';
					}
					else
					{
						push @bsr_atnm,'SC';
					}
					$fflag = 1;
			

		         	}
		
	         	}
	
			if($fpdb =~ /^HETATM/)
			{
		
				if(substr($fpdb,17,3) eq ' CA')
				{
					push @Cax,0.1*sprintf("%f",substr($fpdb,30,8));
					push @Cay,0.1*sprintf("%f",substr($fpdb,38,8));
					push @Caz,0.1*sprintf("%f",substr($fpdb,46,8));
					push @Can,sprintf("%d",substr($fpdb,22,4));
					push @Cachid,substr($fpdb,21,1);
			
				
				}
			
			}
			if($fpdb =~ /^ENDMDL/)
			{
					last;
			}

	
		}
		close FPDB;
	}


	# create .itp file
	my $op_itp_file = $ARGV[0].'.itp';
	open (F1, "> $op_itp_file")or die "cannot open the file $!";


	printf F1 "[moleculetype]\n;molname	exclusions\nProtein  1\n\n\n";


	printf F1 "[atoms]\n ;   nr       type  resnr residue  atom   cgnr     charge\n";
	open ( FH, "< amino_acid_3letter_code.txt") or die "cannot open the file $!";
	open ( FG, "< $atom_type_file") or die "cannot open the file $!";

	$i = 1;
	while(<FH>)
	{
	      my $line = <FG>;;
	      my $str1 = substr($_,0,3);
	      my $str2 = substr($line,0,3);
	      my $i1 = sprintf("%3d",$i);     
		
	   printf F1 "$i1     $str2    $i1      $str1     $str2    $i1      0.000\n";
	      $i = $i + 1;	
	}    
	

	if($metal_flg == 1 && $metal_atom eq'CU')
	{
		my $i11 = sprintf("%3d",$i); 
		my $sizeCu= scalar(@Cun);
		for(my $cc=0; $cc<$sizeCu; $cc++)
		{
			printf F1 "$i11      CU     $i11      CU      CU    $i11      0.000\n";   
			$i11 = $i11+1;
		}

	}
	if($metal_flg == 1 && $metal_atom eq'ZN')
	{
		my $i11 = sprintf("%3d",$i); 
		my $sizeZn= scalar(@Znn);
		for(my $cc=0; $cc<$sizeZn; $cc++)
		{
			printf F1 "$i11      ZN     $i11      ZN      ZN    $i11      0.000\n";   
			$i11 = $i11+1;
		}
	}
	if($metal_flg == 1 && $metal_atom eq'FE')
	{
		my $i11 = sprintf("%3d",$i); 
		my $sizeFe= scalar(@Fen);
		for(my $cc=0; $cc<$sizeFe; $cc++)
		{
			printf F1 "$i11      FE     $i11      FE      FE    $i11      0.000\n";   
			$i11 = $i11+1;
		}
	}
	if($metal_flg == 1 && $metal_atom eq'MG')
	{
		my $i11 = sprintf("%3d",$i); 
		my $sizeMg= scalar(@Mgn);
		for(my $cc=0; $cc<$sizeMg; $cc++)
		{
			printf F1 "$i11      MG     $i11      MG      MG    $i11      0.000\n";   
			$i11 = $i11+1;
		}
	}  
	if($metal_flg == 1 && $metal_atom eq'CA')
	{
		my $i11 = sprintf("%3d",$i); 
		my $sizeCa= scalar(@Can);
		for(my $cc=0; $cc<$sizeCa; $cc++)
		{
			printf F1 "$i11      CA     $i11      CA      CA    $i11      0.000\n";   
			$i11 = $i11+1;
		}
	}
	
	close FH;
	close FG;


	printf F1 "\n\n[ bonds ]\n; i j func-type (here tabulated 8) table number Kb=1 (for tabulated) \n ; virtual bond C-alpha - C-alpha\n"; 


	for( $i = 1;$i<($atom_number);$i++)
	{
		my $i1 = sprintf("%3d",$i);   
	       	my $i2 = sprintf("%3d",$i + 1);  
		my $flag_bond = 0;
		for(my $cct=0; $cct<scalar(@ChEnd); $cct++)
		{
		
			if($i2 == $ChStart[$cct])
			{
						
				$flag_bond = 1; 
			}
		}
		if($flag_bond == 0)
		{
			     
			if($bond_type[$i1-1] == 1)    
			{ printf F1 "$i1    $i2     8       0      1      ;\n"; }
			else
			{  printf F1 "$i1    $i2     8       1      1      ;\n"; }
	
		}	
	

	}

	printf F1 "\n; Disulphide bond\n";

	my @d_nm;
	my @resid;
	my @chain_ty;

	    
	open ( FH, "< $rinfile") or die "cannot open the file $!";
	while(<FH>)
	{
	    
	    my $disul_nm=substr($_,79,1);
	    if ($disul_nm ne " ")
	    {
		  my $res_no = sprintf("%4d",substr($_,0,4));	
		  my $chain = substr($_,8,1);
		    
		  push @d_nm,$disul_nm;
		  push @resid,$res_no;
		  push @chain_ty,$chain;
		 
	     }
	}
	close FH;

	    
	my $len = @resid;

	for(my $i=0;$i<=$len;$i++)
	{
		 
	    for(my $j=$i+1;$j<$len;$j++) 
	    {
		     
		 if($d_nm[$i] eq $d_nm[$j])
		 {
		         
		           printf F1 "$resid[$i]    $resid[$j]     8     900      1\n"   
		 }

		 
	    }

	}    


	my @xcord;
	my @ycord;
	my @zcord;
	$i = 0;
	open ( FH, "< $pdbfile") or die "cannot open the file $!";
	open ( F10, "> t.txt") or die "cannot open the file $!";
	while(<FH>)
	{

		 if (/^ATOM/)
		 {
		         if(substr ($_,13,2) eq 'CA' )
		         {
		                my $x=substr($_,30,8);
		                my $x_val = sprintf("%8.3f",$x);
		                push @xcord,0.1*$x_val;
		                my $y=substr($_,38,8);
		                my $y_val = sprintf("%8.3f",$y);
		                push @ycord,0.1*$y_val;
		                my $z=substr($_,46,8);
		                my $z_val = sprintf("%8.3f",$z); 
		                push @zcord,0.1*$z_val;
		                printf F10 "%8.3f%8.3f%8.3f  %d\n" ,$xcord[$i],$ycord[$i],$zcord[$i],$i; 
		                $i = $i+1;
		               
		                
	 			
		         }


	       }
		if(/^ENDMDL/)
		{
			last;
		}
	}

	
	if($metal_flg == 1 && $metal_atom eq 'CU')
	{
		open ( F23, "> m.txt") or die "cannot open the file $!";
                my $size = scalar(@cu_resno);
                my $size2= scalar(@Cux);


                open (F22, "> metal.txt")or die "cannot open the file $!";
                my $cu_cgmm_res_no = $totalnores;
                my $fl = 0;
                for(my $cn=0; $cn<$size2; $cn++)
                {
                        print F23 "$Cun[$cn]  $Cux[$cn]  $Cuy[$cn]  $Cuz[$cn]\n";
                        $cu_cgmm_res_no = $cu_cgmm_res_no + 1;
                        for(my $cntt=0; $cntt<$size; $cntt++)
                        {
                                my $Sz= scalar(@ori_resid_no);
                                for(my $cntt2=0; $cntt2<$Sz; $cntt2++)
                                {
                                       
                                        if($ori_resid_no[$cntt2] == $bsr_resno[$cntt] && $bsr_chid[$cntt] eq $ori_chid[$cntt2])
                                        {
                                                my $str = $cu_chid[$cntt]. $bsr_chid[$cntt];
						my $str1 = $Cuchid[$cn]. $bsr_chid[$cntt];

						if($Cun[$cn] == $cu_resno[$cntt]&& $str eq $str1)
                                                {
                                                        
                                                       printf F22 "%4d%8.3f%8.3f%8.3f%4d%5s%8.3f%8.3f%8.3f\n",$cu_cgmm_res_no,$Cux[$cn],$Cuy[$cn],$Cuz[$cn],$cgmm_resid_no[$cntt2], $bsr_resnm[$cntt],$xcord[$cgmm_resid_no[$cntt2]-1],$ycord[$cgmm_resid_no[$cntt2]-1],$zcord[$cgmm_resid_no[$cntt2]-1];
                                                        
                                                }
                                        }
                                }
                        

                        }
                }
                close F22;
                close F23;

		system("awk '!a[\$0]++' ./metal.txt > ./Mtemp.txt");
		system("cp ./Mtemp.txt ./metal.txt");
		system("rm ./Mtemp.txt");

                printf F1 "\n; metal bond\n";
                open( FM, "< metal.txt") or die "cannot open the file $!";
                while(my $fmetal = <FM>)
                {
                        my @metal_line = split(' ',$fmetal);
                        
                        my $dxm = $metal_line[1] - $metal_line[6];  
                        my $dym = $metal_line[2] - $metal_line[7];  
                        my $dzm = $metal_line[3] - $metal_line[8]; 
                        
                        my $dist_metal = sqrt($dxm * $dxm + $dym * $dym + $dzm * $dzm);

                        if($dist_metal < 6)
                        { 
                                printf F1 "$metal_line[0]   $metal_line[4]     8       903      1      ;\n";
                        }
                        else
                        {
                                printf F1 "$metal_line[0]   $metal_line[4]     8       902      1      ;\n";

                        }
                }
                close FM;
		

            }
		
	    
	    if($metal_flg == 1 && $metal_atom eq 'ZN')
	    {
		open ( F23, "> m.txt") or die "cannot open the file $!";
		my $size = scalar(@zn_resno);
		my $size2= scalar(@Znx);
		

		open (F22, "> metal.txt")or die "cannot open the file $!";
		my $zn_cgmm_res_no = $totalnores;
		
		for(my $cn=0; $cn<$size2; $cn++)
		{
			
			print F23 "$Znn[$cn]  $Znx[$cn]  $Zny[$cn]  $Znz[$cn]\n";
			$zn_cgmm_res_no = $zn_cgmm_res_no + 1;
			
			for(my $cntt=0; $cntt<$size; $cntt++)
			{
				
				my $Sz= scalar(@ori_resid_no);
				for(my $cntt2=0; $cntt2<$Sz; $cntt2++)
				{
				        
					if($ori_resid_no[$cntt2] == $bsr_resno[$cntt]&& $bsr_chid[$cntt] eq $ori_chid[$cntt2])
					{
						my $str = $zn_chid[$cntt]. $bsr_chid[$cntt];
						my $str1 = $Znchid[$cn]. $bsr_chid[$cntt];
						
						if($Znn[$cn] == $zn_resno[$cntt] && $str eq $str1)
						{
							
								
								printf F22 "%4d%8.3f%8.3f%8.3f%4d%5s%8.3f%8.3f%8.3f\n",$zn_cgmm_res_no,$Znx[$cn],$Zny[$cn],$Znz[$cn],$cgmm_resid_no[$cntt2], $bsr_resnm[$cntt],$xcord[$cgmm_resid_no[$cntt2]-1],$ycord[$cgmm_resid_no[$cntt2]-1],$zcord[$cgmm_resid_no[$cntt2]-1];
							}
						
					}
				}
	

			}
		}
		close F22;
		close F23;
		
		system("awk '!a[\$0]++' ./metal.txt > ./Mtemp.txt");
		system("cp ./Mtemp.txt ./metal.txt");
		system("rm ./Mtemp.txt");
		

		printf F1 "\n; metal bond\n";
		open( FM, "< metal.txt") or die "cannot open the file $!";
		while(my $fmetal = <FM>)
		{
			my @metal_line = split(' ',$fmetal);
	
			my $dxm = $metal_line[1] - $metal_line[6];  
		    	my $dym = $metal_line[2] - $metal_line[7];  
		    	my $dzm = $metal_line[3] - $metal_line[8]; 
	
			my $dist_metal = sqrt($dxm * $dxm + $dym * $dym + $dzm * $dzm);

			if($dist_metal < 6)
			{ 
				printf F1 "$metal_line[0]   $metal_line[4]     8       901      1      ;\n";
			}
			else
			{
				printf F1 "$metal_line[0]   $metal_line[4]     8       902      1      ;\n";

			}
		}
		close FM;

	}
	
	if($metal_flg == 1 && $metal_atom eq 'FE')
	{
		
		open ( F23, "> m.txt") or die "cannot open the file $!";
                my $size = scalar(@fe_resno);
                my $size2= scalar(@Fex);
		
		

                open (F22, "> metal.txt")or die "cannot open the file $!";
                my $fe_cgmm_res_no = $totalnores;
               
                for(my $cn=0; $cn<$size2; $cn++)
                {
                        
			print F23 "$Fen[$cn]  $Fex[$cn]  $Fey[$cn]  $Fez[$cn]\n";
                        $fe_cgmm_res_no = $fe_cgmm_res_no + 1;
                        for(my $cntt=0; $cntt<$size; $cntt++)
                        {
				
                                my $Sz= scalar(@ori_resid_no);
                                for(my $cntt2=0; $cntt2<$Sz; $cntt2++)
                                {
                                                                
                                        if($ori_resid_no[$cntt2] == $bsr_resno[$cntt] && $bsr_chid[$cntt] eq $ori_chid[$cntt2])
                                        {
                                                my $str = $fe_chid[$cntt]. $bsr_chid[$cntt];
						my $str1 = $Fechid[$cn]. $bsr_chid[$cntt];
						

						if($Fen[$cn] == $fe_resno[$cntt]&& $str eq $str1)
                                                {
                                                       
                                                                printf F22 "%4d%8.3f%8.3f%8.3f%4d%5s%8.3f%8.3f%8.3f\n",$fe_cgmm_res_no,$Fex[$cn],$Fey[$cn],$Fez[$cn],$cgmm_resid_no[$cntt2], $bsr_resnm[$cntt],$xcord[$cgmm_resid_no[$cntt2]-1],$ycord[$cgmm_resid_no[$cntt2]-1],$zcord[$cgmm_resid_no[$cntt2]-1];
                                                       
                                                }
                                        }
                                }
                        

                        }
                }
                close F22;
                close F23;
		
		system("awk '!a[\$0]++' ./metal.txt > ./Mtemp.txt");
		system("cp ./Mtemp.txt ./metal.txt");
		system("rm ./Mtemp.txt");
	
                printf F1 "\n; metal bond\n";
                open( FM, "< metal.txt") or die "cannot open the file $!";
                while(my $fmetal = <FM>)
                {
                        my @metal_line = split(' ',$fmetal);
                        
                        my $dxm = $metal_line[1] - $metal_line[6];  
                        my $dym = $metal_line[2] - $metal_line[7];  
                        my $dzm = $metal_line[3] - $metal_line[8]; 
                        
                        my $dist_metal = sqrt($dxm * $dxm + $dym * $dym + $dzm * $dzm);

                        printf F1 "$metal_line[0]   $metal_line[4]     8       906      1      ;\n";
                        
                }
                close FM;


	}

	if($metal_flg == 1 && $metal_atom eq'MG')
	{
		open ( F23, "> m.txt") or die "cannot open the file $!";
		my $size = scalar(@mg_resno);
		my $size2= scalar(@Mgx);


		open (F22, "> metal.txt")or die "cannot open the file $!";
		my $mg_cgmm_res_no = $totalnores;
		


		for(my $cn=0; $cn<$size2; $cn++)
		{
			print F23 "$Mgn[$cn]  $Mgx[$cn]  $Mgy[$cn]  $Mgz[$cn]\n";
			$mg_cgmm_res_no = $mg_cgmm_res_no + 1;


			for(my $cntt=0; $cntt<$size; $cntt++)
			{
				my $Sz= scalar(@ori_resid_no);
				for(my $cntt2=0; $cntt2<$Sz; $cntt2++)
				{
						
					if($ori_resid_no[$cntt2] == $bsr_resno[$cntt]&& $bsr_chid[$cntt] eq $ori_chid[$cntt2])
					{
						
						my $str = $mg_chid[$cntt]. $bsr_chid[$cntt];
						my $str1 = $Mgchid[$cn]. $bsr_chid[$cntt];
								
						if($Mgn[$cn] == $mg_resno[$cntt]&& $str eq $str1)
						{
							
						
									printf F22 "%4d%8.3f%8.3f%8.3f%4d%5s%8.3f%8.3f%8.3f%5s\n",$mg_cgmm_res_no,$Mgx[$cn],$Mgy[$cn],$Mgz[$cn],$cgmm_resid_no[$cntt2], $bsr_resnm[$cntt],$xcord[$cgmm_resid_no[$cntt2]-1],$ycord[$cgmm_resid_no[$cntt2]-1],$zcord[$cgmm_resid_no[$cntt2]-1],$bsr_atnm[$cntt]; 					
						
						
						
							
						}
						
					}
				}
	

			}
		}
		close F22;
		close F23;
	
		system("awk '!a[\$0]++' ./metal.txt > ./Mtemp.txt");
		system("cp ./Mtemp.txt ./metal.txt");
		system("rm ./Mtemp.txt");
		

		printf F1 "\n; metal bond\n";
		open( FM, "< metal.txt") or die "cannot open the file $!";
		while(my $fmetal = <FM>)
		{
			my @metal_line = split(' ',$fmetal);
	
			my $dxm = $metal_line[1] - $metal_line[6];  
		    	my $dym = $metal_line[2] - $metal_line[7];  
		    	my $dzm = $metal_line[3] - $metal_line[8]; 
	
			my $dist_metal = sqrt($dxm * $dxm + $dym * $dym + $dzm * $dzm);

	
			printf F1 "$metal_line[0]   $metal_line[4]     8       905      1      ;\n";
	
	
		}
		close FM;

	}
	
	if($metal_flg == 1 && $metal_atom eq 'CA')	 
	{
			
		open ( F23, "> m.txt") or die "cannot open the file $!";
		my $size = scalar(@ca_resno);
		my $size2= scalar(@Cax);
		

		open (F22, "> metal.txt")or die "cannot open the file $!";
		my $ca_cgmm_res_no = $totalnores;
		
		
		for(my $cn=0; $cn<$size2; $cn++)
		{
			print F23 "$Can[$cn]  $Cax[$cn]  $Cay[$cn]  $Caz[$cn]\n";
			$ca_cgmm_res_no = $ca_cgmm_res_no + 1;

			

			for(my $cntt=0; $cntt<$size; $cntt++)
			{
				my $Sz= scalar(@ori_resid_no);
				for(my $cntt2=0; $cntt2<$Sz; $cntt2++)
				{
						
					if($ori_resid_no[$cntt2] == $bsr_resno[$cntt] && $bsr_chid[$cntt] eq $ori_chid[$cntt2])
					{
						my $str = $ca_chid[$cntt]. $bsr_chid[$cntt];
						my $str1 = $Cachid[$cn]. $bsr_chid[$cntt];
		
						if($Can[$cn] == $ca_resno[$cntt] && $str eq $str1)
						{
							
							printf F22 "%4d%8.3f%8.3f%8.3f%4d%5s%8.3f%8.3f%8.3f%5s\n",$ca_cgmm_res_no,$Cax[$cn],$Cay[$cn],$Caz[$cn],$cgmm_resid_no[$cntt2], $bsr_resnm[$cntt],$xcord[$cgmm_resid_no[$cntt2]-1],$ycord[$cgmm_resid_no[$cntt2]-1],$zcord[$cgmm_resid_no[$cntt2]-1],$bsr_atnm[$cntt]; 					
										
							
						}
						
					}
				}
	

			}
		}
		close F22;
		close F23;

		system("awk '!a[\$0]++' ./metal.txt > ./Mtemp.txt");
		system("cp ./Mtemp.txt ./metal.txt");
		system("rm ./Mtemp.txt");

		printf F1 "\n; metal bond\n";
		open( FM, "< metal.txt") or die "cannot open the file $!";
		while(my $fmetal = <FM>)
		{
			my @metal_line = split(' ',$fmetal);
	
			my $dxm = $metal_line[1] - $metal_line[6];  
    			my $dym = $metal_line[2] - $metal_line[7];  
    			my $dzm = $metal_line[3] - $metal_line[8]; 
	
			my $dist_metal = sqrt($dxm * $dxm + $dym * $dym + $dzm * $dzm);

	
			printf F1 "$metal_line[0]   $metal_line[4]     8       904      1      ;\n";
	
	
		}
		close FM;
	}
	
        
	close FH;
	close F10;
	printf F1 "\n; nonbonded\n";
	open( FG, "< nonbonded_pair.txt") or die "cannot open the file $!";

	
	while(<FG>)
	{
	    
	     my $nb1 = substr($_,0,5);
	     my $nb2 = substr($_,5,5);
	     my $nb3 = substr($_,11,2);
	     my $nb3_prime = substr($_,12,1).substr($_,11,1);

	     my $num_nb1 = sprintf("%6d",$nb1)-1;
	     my $num_nb2 = sprintf("%6d",$nb2)-1;
	  
	 
	    my $dx = $xcord[$num_nb1] - $xcord[$num_nb2];  
	    my $dy = $ycord[$num_nb1] - $ycord[$num_nb2];  
	    my $dz = $zcord[$num_nb1] - $zcord[$num_nb2];  
	    
	  
		 
	    
	    

	    my $distance = sqrt($dx * $dx + $dy * $dy + $dz * $dz);
	    
	    open( FH, "< nbclust_numb.txt") or die "cannot open the file $!"; 
	   
	    
	     while(<FH>)
	     {
		   
	    
		    my $bd1 = 0.1*5.29;
		    my $bd2 = 0.1*5.89;
		    my $bd3 = 0.1*7.7;
		  
		  my $nb4 = substr($_,0,2);  
		  if($nb3 eq $nb4 || $nb3_prime eq $nb4)
		  {
		       my $nb5 = sprintf("%5d",substr($_,3,3));

		        
	 
		      if( $distance < $bd3)
		       {
		       	       my $nb5_1 = $nb5 + 300;	
		               printf F1 "$nb1   $nb2     8       $nb5_1      1      ;\n";     
		       }
		       else
		       { 
			       if($distance < 1.2)
		               {
			
		                   printf F1 "$nb1   $nb2     8     $nb5      1      ;\n"; 
		               }
		       }		
		  }
	     }
	     close FH;
	     
	}


	close FG;


	reset 'xcord';
	reset 'ycord';
	reset 'zcord';


	close FG;



	printf F1 "\n\n[ angles ]\n; i   j   k   func-type(here tabulated 8)   tablenumber   constant=1(for tbulated only) \n ; virtual bond angle C-alpha - C-alpha - C-alpha\n"; 
	open( FG, "< triplet.txt") or die "cannot open the file $!";

	while(<FG>)
	{
	     my $a1 = substr($_,0,4);
	     my $a2 = substr($_,5,4);
	     my $a3 = substr($_,10,4);
	     my $a4 = substr($_,15,3);

	    
	     open( FH, "< triplet2cluster.txt") or die "cannot open the file $!"; 
	     while(<FH>)
	     {
		  my $a5 = substr($_,0,3);  
		  
		  if($a4 eq $a5 )
		  {
		      
		       my $a6 = substr($_,4,3);
		       my $a1_no = sprintf("%d",$a1);
		       my $a2_no = sprintf("%d",$a2);
		       my $a3_no = sprintf("%d",$a3);
		       my $flag_ang = 0;
		       for(my $cct=0; $cct<scalar(@ChEnd); $cct++)
		       {
			    if(($a2_no == $ChEnd[$cct] && $a3_no == $ChStart[$cct]) || ($a1_no == $ChEnd[$cct] && $a2_no == $ChStart[$cct]))
			    {
				$flag_ang = 1; 
			    }
		       }
		       for(my $cct1 = 0; $cct1<scalar(@domainStart); $cct1++)
		       {
			    
			    if($a3_no >= $domainStart[$cct1] && $a1_no <= $domainEnd[$cct1])
			    {
				$flag_ang = 1;
			    }

		       }		
		       if($flag_ang == 0)             
		       {printf F1 "$a1   $a2    $a3        8     $a6      1    ;\n";     }
		  
		  }
	     }
	     close FH;
	}


	close FG;

	my $metal_size;
	my @metal_data;

	
	if($metal_flg == 1 && $metal_atom eq 'CU') 
	{
		open( FM, "< metal.txt") or die "cannot open the file $!";
		@metal_data=<FM>;
		close FM;
		$metal_size = scalar(@metal_data);
		for(my $mcn=0; $mcn<$metal_size; $mcn++)
		{
			my @clm1 = split(' ',$metal_data[$mcn]);
			for(my $mcn1=$mcn+1; $mcn1<$metal_size; $mcn1++)
			{
				my @clm2 = split(' ',$metal_data[$mcn1]);
				my $aapair = $clm1[5].$clm2[5];
		
				if($clm1[0] == $clm2[0])
				{
					my $tab_no=990;
			
					printf F1 "%3d%6d%7d        8     $tab_no      1    ;\n", $clm1[4],$clm1[0],$clm2[4];    
				}
			}

		}
	}

	if($metal_flg == 1 && $metal_atom eq 'ZN') 
	{
		open( FM, "< metal.txt") or die "cannot open the file $!";
		@metal_data=<FM>;
		close FM;
		$metal_size = scalar(@metal_data);
		for(my $mcn=0; $mcn<$metal_size; $mcn++)
		{
			my @clm1 = split(' ',$metal_data[$mcn]);
			for(my $mcn1=$mcn+1; $mcn1<$metal_size; $mcn1++)
			{
				my @clm2 = split(' ',$metal_data[$mcn1]);
				my $aapair = $clm1[5].$clm2[5];
		
				if($clm1[0] == $clm2[0])
				{
					my $tab_no=0;
					if($aapair eq 'CYSCYS' || $aapair eq 'CYSHIS' || $aapair eq 'HISHIS')
					{
						if($aapair eq 'CYSCYS')
						{		
							$tab_no = 999;
						}
						if($aapair eq 'CYSHIS')
						{
							$tab_no = 998;
						}
						if($aapair eq 'CYSHIS')
						{
							$tab_no = 995;
						}

					}
					else
					{
						$tab_no = 991;
					}
			
					printf F1 "%3d%6d%7d        8     $tab_no      1    ;\n", $clm1[4],$clm1[0],$clm2[4];    
				}
			}

		}
	}
	if($metal_flg == 1 && $metal_atom eq 'FE') 
	{
		open( FM, "< metal.txt") or die "cannot open the file $!";
		@metal_data=<FM>;
		close FM;
		$metal_size = scalar(@metal_data);
		for(my $mcn=0; $mcn<$metal_size; $mcn++)
		{
			my @clm1 = split(' ',$metal_data[$mcn]);
			for(my $mcn1=$mcn+1; $mcn1<$metal_size; $mcn1++)
			{
				my @clm2 = split(' ',$metal_data[$mcn1]);
				my $aapair = $clm1[5].$clm2[5];
		
				if($clm1[0] == $clm2[0])
				{
					my $tab_no= 889;
			
					printf F1 "%3d%6d%7d        8     $tab_no      1    ;\n", $clm1[4],$clm1[0],$clm2[4];    
				}
			}

		}

	}
	if($metal_flg == 1 && $metal_atom eq 'MG')
	{
		open( FM, "< metal.txt") or die "cannot open the file $!";
		@metal_data=<FM>;
		close FM;

		$metal_size = scalar(@metal_data);
		for(my $mcn=0; $mcn<$metal_size; $mcn++)
		{
			my @clm1 = split(' ',$metal_data[$mcn]);
			for(my $mcn1=$mcn+1; $mcn1<$metal_size; $mcn1++)
			{
				my @clm2 = split(' ',$metal_data[$mcn1]);
				my $aapair = $clm1[5].$clm2[5];
				my $str_ms = $clm1[9].$clm2[9];
		
				my $tab_no=0;
				if($clm1[0] == $clm2[0])
				{
					if($str_ms eq 'MCMC' || $str_ms eq 'MCSC' || $str_ms eq 'SCMC'|| $str_ms eq 'SCSC')
					{
						if($str_ms eq 'MCMC')
						{
							$tab_no = 982;
						}
						if($str_ms eq 'MCSC'|| $str_ms eq 'SCMC')
						{
							$tab_no = 981;
						}
						if($str_ms eq 'SCSC')
						{
							$tab_no = 980;
						}
					}
		
				printf F1 "%3d%6d%7d        8     $tab_no      1    ;\n", $clm1[4],$clm1[0],$clm2[4];    
				}
			}

		}


	}
	if($metal_flg == 1 && $metal_atom eq 'CA') 
	{
		open( FM, "< metal.txt") or die "cannot open the file $!";
		@metal_data=<FM>;
		close FM;
		my $flag=0;
		$metal_size = scalar(@metal_data);
		for(my $mcn=0; $mcn<$metal_size; $mcn++)
		{
			my @clm1 = split(' ',$metal_data[$mcn]);
			for(my $mcn1=$mcn+1; $mcn1<$metal_size; $mcn1++)
			{
				my @clm2 = split(' ',$metal_data[$mcn1]);
				my $aapair = $clm1[5].$clm2[5];
				my $str_ms = $clm1[9].$clm2[9];
		
				my $tab_no=0;
				if($clm1[0] == $clm2[0])
				{
					if($aapair eq 'ASNASP' ||  $aapair eq 'ASPASN' || $aapair eq 'ASPASP' || $aapair eq 'ASPGLU' || $aapair eq 'GLUASP' || $aapair eq 'GLUGLU')
					{				
						if($aapair eq 'ASNASP' ||  $aapair eq 'ASPASN')
						{
							$tab_no = 988; $flag = 1;
						}
						if($aapair eq 'ASPASP')
						{
							$tab_no = 987; $flag = 1;
						}
						if($aapair eq 'ASPGLU')
						{
							$tab_no = 986; $flag = 1;

						}
						if($aapair eq 'GLUASP')
						{
							$tab_no = 986; $flag = 1;
						}
						if($aapair eq 'GLUGLU')
						{
							$tab_no = 985; $flag = 1;

						}
					}
					else
					{			
						if($flag == 0)
						{
							if($str_ms eq 'MCMC' || $str_ms eq 'MCSC' || $str_ms eq 'SCMC'|| $str_ms eq 'SCSC')
							{
								if($str_ms eq 'MCMC')
								{
									$tab_no = 985;
								}
								if($str_ms eq 'MCSC'|| $str_ms eq 'SCMC')
								{
									$tab_no = 984;
								}
								if($str_ms eq 'SCSC')
								{
									$tab_no = 983;
								}

							}

						}
						else
						{
							$tab_no = 985;

						}
					}
		
					printf F1 "%3d%6d%7d        8     $tab_no      1    ;\n", $clm1[4],$clm1[0],$clm2[4];    
				}
			}

		    }	


	}
	

	
	
	printf F1 "\n\n[ dihedrals ]\n; i   j   k   n   func-type(here tabulated 8)   tablenumber   constant=1(for tbulated only) \n ; virtual bond angle C-alpha - C-alpha - C-alpha - C-alpha\n"; 
	open( FG, "< quartet.txt") or die "cannot open the file $!";

	while(my $line1=<FG>)
	{
	     my $d1 = substr($line1,0,4);
	     my $d2 = substr($line1,5,4);
	     my $d3 = substr($line1,10,4);
	     my $d4 = substr($line1,15,4);
	     my $d5 = substr($line1,20,4);
	     
	     my $d1_no = sprintf("%d",$d1) - 1;
	     my $d2_no = sprintf("%d",$d2) - 1;
	     my $d3_no = sprintf("%d",$d3) - 1;
	     my $d4_no = sprintf("%d",$d4) - 1;

	     my $D1_no = sprintf("%d",$d1) ;
	     my $D2_no = sprintf("%d",$d2);
	     my $D3_no = sprintf("%d",$d3);
	     my $D4_no = sprintf("%d",$d4);
	
	     my $flag_dih = 0;	
	     
	     if($atom[$d1_no] eq 'G' | $atom[$d2_no] eq 'G' | $atom[$d3_no] eq 'G' | $atom[$d4_no] eq 'G' )
	     {
		    open( FH, "<quartet2clusterGLY.txt") or die "cannot open the file $!"; 
	     	    while(my $line2=<FH>)
	     	    {
		  	my $d6 = substr($line2,0,4);  
		  	if($d5 eq $d6 )
		  	{
		       
		       		my $d7 = sprintf("%d",substr($line2,5,3)) + 400;
				for(my $cct=0; $cct<scalar(@ChEnd); $cct++)
		       		{	
			    		if(($D1_no == $ChEnd[$cct] && $D2_no == $ChStart[$cct]) || ($D2_no ==  $ChEnd[$cct] && $D3_no == $ChStart[$cct]) || ($D3_no ==  $ChEnd[$cct] && $D4_no == $ChStart[$cct]))
			    		{
						$flag_dih = 1; 
			    		}
		       		}
				for(my $cct1 = 0; $cct1<scalar(@domainStart); $cct1++)
		       		{
			    
			    		if($D4_no >= $domainStart[$cct1] && $D1_no <= $domainEnd[$cct1])
			    		{
						$flag_dih = 1;
			    		}

		       		}	
		                if($flag_dih == 0)
				{   
		       		printf F1 "$d1    $d2    $d3    $d4      8     $d7      1      ;\n";  
		       		reset $d7;
				}
		  
		  	}
		  	reset $d6;
	     	    }
	     	    close FH;

	     }
	     else
	     {		
	     	    open( FH, "<quartet2cluster.txt") or die "cannot open the file $!"; 
	     	    while(my $line2=<FH>)
	     	    {
		  	my $d6 = substr($line2,0,4);  
		  	if($d5 eq $d6 )
		  	{
		       
		       		my $d7 = substr($line2,5,3);
		                for(my $cct=0; $cct<scalar(@ChEnd); $cct++)
		       		{	
			    		if(($D1_no == $ChEnd[$cct] && $D2_no == $ChStart[$cct]) || ($D2_no ==  $ChEnd[$cct] && $D3_no == $ChStart[$cct]) || ($D3_no ==  $ChEnd[$cct] && $D4_no == $ChStart[$cct]))
			    		{
						$flag_dih = 1; 
			    		}
		       		}
				for(my $cct1 = 0; $cct1<scalar(@domainStart); $cct1++)
		       		{
			    
			    		if($D4_no >= $domainStart[$cct1] && $D1_no <= $domainEnd[$cct1])
			    		{
						$flag_dih = 1;
			    		}

		       		}	
			
		                if($flag_dih == 0)
				{      
		       			printf F1 "$d1    $d2    $d3    $d4      8     $d7      1      ;\n";     
		       			reset $d7;
				}
		  
		  	}
		  	reset $d6;
	     	    }
	     	    close FH;
	     }
	     
	}


	close FG;

	

	if($metal_flg == 1  && $metal_atom eq 'CU') 
	{
		for(my $mcn=0; $mcn<$metal_size; $mcn++)
		{
			my @clm1 = split(' ',$metal_data[$mcn]);
			for(my $mcn1=0; $mcn1<$metal_size; $mcn1++)
			{
				my @clm2 = split(' ',$metal_data[$mcn1]);		
				for(my $mcn2 = $mcn1+1; $mcn2<$metal_size; $mcn2++)
				{
					my @clm3 = split(' ',$metal_data[$mcn2]);
					if($mcn != $mcn1 && $mcn != $mcn2)
					{
						my $aatrpl = $clm1[5].$clm2[5].$clm3[5];
				
						if($clm1[0] == $clm2[0] && $clm3[0] == $clm2[0] )
						{				
							my $tab_no=989;
					
							printf F1 "%3d%7d%7d%7d      8     $tab_no      1    ;\n", $clm2[4],$clm1[0],$clm1[4],$clm3[4]; 
						}
					}
				}   
			}

		}

	}
	if($metal_flg == 1  && $metal_atom eq 'ZN') 
	{
		for(my $mcn=0; $mcn<$metal_size; $mcn++)
		{
			my @clm1 = split(' ',$metal_data[$mcn]);
			for(my $mcn1=0; $mcn1<$metal_size; $mcn1++)
			{
				my @clm2 = split(' ',$metal_data[$mcn1]);		
				for(my $mcn2 = $mcn1+1; $mcn2<$metal_size; $mcn2++)
				{
					my @clm3 = split(' ',$metal_data[$mcn2]);
					if($mcn != $mcn1 && $mcn != $mcn2)
					{
						my $aatrpl = $clm1[5].$clm2[5].$clm3[5];
				
						if($clm1[0] == $clm2[0] && $clm3[0] == $clm2[0] )
						{				
							my $tab_no=0;
							if($aatrpl eq 'CYSCYSCYS' || $aatrpl eq 'HISHISHIS')
							{
								if($aatrpl eq 'CYSCYSCYS')
								{
									$tab_no = 999;
								}
								if($aatrpl eq 'HISHISHIS')
								{
									$tab_no = 995;
								}
							}
							else
							{
								$tab_no = 990;
							}
					
						printf F1 "%3d%7d%7d%7d      8     $tab_no      1    ;\n", $clm2[4],$clm1[0],$clm1[4],$clm3[4]; 
						}
					}
				}   
			}

		}



	}
	if($metal_flg == 1 && $metal_atom eq 'FE') 
	{
		for(my $mcn=0; $mcn<$metal_size; $mcn++)
		{
			my @clm1 = split(' ',$metal_data[$mcn]);
			for(my $mcn1=0; $mcn1<$metal_size; $mcn1++)
			{
				my @clm2 = split(' ',$metal_data[$mcn1]);		
				for(my $mcn2 = $mcn1+1; $mcn2<$metal_size; $mcn2++)
				{
					my @clm3 = split(' ',$metal_data[$mcn2]);
					if($mcn != $mcn1 && $mcn != $mcn2)
					{
						my $aatrpl = $clm1[5].$clm2[5].$clm3[5];
				
						if($clm1[0] == $clm2[0] && $clm3[0] == $clm2[0] )
						{				
							my $tab_no=888;
					
							printf F1 "%3d%7d%7d%7d      8     $tab_no      1    ;\n", $clm2[4],$clm1[0],$clm1[4],$clm3[4]; 
						}
					}
				}   
			}

		}


	}
	
	if($metal_flg == 1 && $metal_atom eq 'MG')
	{
		for(my $mcn=0; $mcn<$metal_size; $mcn++)
		{
			my @clm1 = split(' ',$metal_data[$mcn]);
			for(my $mcn1=0; $mcn1<$metal_size; $mcn1++)
			{
				my @clm2 = split(' ',$metal_data[$mcn1]);		
				for(my $mcn2 = $mcn1+1; $mcn2<$metal_size; $mcn2++)
				{
					my @clm3 = split(' ',$metal_data[$mcn2]);
					if($mcn != $mcn1 && $mcn != $mcn2)
					{
						my $aatrpl = $clm1[5].$clm2[5].$clm3[5];
						my $str_scc = $clm1[9].$clm2[9].$clm2[9];
						if($clm1[0] == $clm2[0] && $clm3[0] == $clm2[0] )
						{
				
							my $tab_no=0;
							if($str_scc eq 'SCSCMC' || $str_scc eq 'MCSCSC' )
							{
								$tab_no = 984; 
							}
							if($str_scc eq 'SCMCSC' || $str_scc eq 'SCSCMC' )
							{
								$tab_no = 983; 
					
							}
							if($str_scc eq 'MCMCMC' || $str_scc eq 'MCMCSC' || $str_scc eq 'MCSCMC' || $str_scc eq 'SCMCMC')
							{
								$tab_no = 982; 
							}
							if($tab_no != 0)
							{
								printf F1 "%3d%7d%7d%7d      8     $tab_no      1    ;\n", $clm2[4],$clm1[0],$clm1[4],$clm3[4]; 						}	
						}
					}
				}   
			}

		}

	}

	if($metal_flg == 1  && $metal_atom eq 'CA') 
	{

		for(my $mcn=0; $mcn<$metal_size; $mcn++)
		{
			my @clm1 = split(' ',$metal_data[$mcn]);
			for(my $mcn1=0; $mcn1<$metal_size; $mcn1++)
			{
				my @clm2 = split(' ',$metal_data[$mcn1]);		
				for(my $mcn2 = $mcn1+1; $mcn2<$metal_size; $mcn2++)
				{
					my @clm3 = split(' ',$metal_data[$mcn2]);
					if($mcn != $mcn1 && $mcn != $mcn2)
					{
						my $aatrpl = $clm1[5].$clm2[5].$clm3[5];
						my $str_scc = $clm1[9].$clm2[9].$clm2[9];
						if($clm1[0] == $clm2[0] && $clm3[0] == $clm2[0] )
						{
				
							my $tab_no=0;
							if($str_scc eq 'SCSCMC' || $str_scc eq 'MCSCSC' || $str_scc eq 'SCMCSC')
							{
								$tab_no = 988; 
							}
							if($str_scc eq 'SCSCSC' || $str_scc eq 'SCMCMC' )
							{
								$tab_no = 987; 
					
							}
							if($str_scc eq 'MCSCMC' || $str_scc eq 'MCMCSC')
							{
								$tab_no = 986; 
							}
							if($str_scc eq 'MCMCMC' )
							{
								$tab_no = 985; 
							}
							printf F1 "%3d%7d%7d%7d      8     $tab_no      1    ;\n", $clm2[4],$clm1[0],$clm1[4],$clm3[4]; 
						}
					}
				}   
			}

		}
	}	
	

	close F1;

	# create .top file
	my $op_top_file = $ARGV[0].'.top';
	open (F2, "> $op_top_file")or die "cannot open the file $!";
	my $ff = '"' . 'cgmm.itp' .'"';
	my $itp = '"'. $op_itp_file . '"';
	my $ucpdb = uc($ARGV[0]);
	printf F2 "#include $ff\n#include $itp\n\n[ system ]\n; name\n$ucpdb\n\n\n[ molecules ]\n; name  number\nProtein	1";
	close F2;

	#copy itp and top in a particuler folder
	system("rm *.rin");
	system("rm *.ssd");
	system("rm *.seq");
	system("rm nonbonded_pair.txt quartet.txt t.txt triplet.txt amino_acid_3letter_code.txt c.txt ");


	#create .gro file
	# create .gro file
	
	if($metal_flg == 1)
	{
				
		$filename = $ARGV[0]; 
		my $proteinname = $ARGV[0];
		$pdbfile = $ARGV[0]. '.pdb';
		my $grofile = $ARGV[0]. '.gro';
		open ( FH, "< $pdbfile") or die "cannot open the file $!"; 


		my @xcord1;
		my @ycord1;
		my @zcord1;
		my @res_name;
		$i = 0;
		while(<FH>)
		{

		 	if (/^ATOM/)
		 	{
		         	if(substr ($_,13,2) eq 'CA' )
		         	{
		                	my $x=substr($_,30,8);
		                	my $x_val = sprintf("%8.3f",$x);
		                	push @xcord1,0.1*$x_val;
		                	my $y=substr($_,38,8);
		                	my $y_val = sprintf("%8.3f",$y);
		                	push @ycord1,0.1*$y_val;
		                	my $z=substr($_,46,8);
		                	my $z_val = sprintf("%8.3f",$z); 
		                	push @zcord1,0.1*$z_val;
		                	$i = $i+1;
		               
		                	my $res = substr($_,17,3);
		                	push @res_name,$res;
	 			
		         	}


	       		}
		}

		reset 'xcord1';
		reset 'ycord1';
		reset 'zcord1';
		close FH;
		$atom_type_file = ($ARGV[0]).'_atom_type.txt';
		open ( FG, "< $atom_type_file") or die "cannot open the file $!";
		my @at;
		while(<FG>)
		{
			my $a = substr($_,0,3);
			push @at,$a;

		}
		close FG;

		$len = scalar(@at);

		open ( FM, "< m.txt") or die "cannot open the file $!";
		my @metal=<FM>;
		close FM;
		
		my $lenm = scalar(@metal);
		my $len_total = $len + $lenm;

		open (F1, "> $grofile")or die "cannot open the file $!";
		print F1 "$proteinname\n  $len_total\n";
		for($i = 0; $i<$len ; $i+=1)
		{
			printf F1 "%5d%5s%5s%5d%8.3f%8.3f%8.3f\n" , $i+1,$res_name[$i],$at[$i],$i+1,$xcord1[$i],$ycord1[$i],$zcord1[$i];
		}
		my $mion = $metal_atom;
		for(my $ii = 0; $ii < $lenm; $ii++)
		{
			my @term = split(' ',$metal[$ii]);
	
			printf F1 "%5d%5s%5s%5d%8.3f%8.3f%8.3f\n" , $i+1,$mion,$mion,$i+1,$term[1],$term[2],$term[3];
			$i = $i+1;
		}	
		print F1 "  0.00000   0.00000   0.00000";

		close F1;
	}
	else
	{
		$filename = $ARGV[0]; 
		my $proteinname = $ARGV[0];
		$pdbfile = $ARGV[0]. '.pdb';
		my $grofile = $ARGV[0]. '.gro';
		open ( FH, "< $pdbfile") or die "cannot open the file $!"; 


		my @xcord1;
		my @ycord1;
		my @zcord1;
		my @res_name;
		$i = 0;
		while(<FH>)
		{

		 	if (/^ATOM/)
		 	{
		         	if(substr ($_,13,2) eq 'CA' )
		         	{
		                	my $x=substr($_,30,8);
		                	my $x_val = sprintf("%8.3f",$x);
		                	push @xcord1,0.1*$x_val;
		                	my $y=substr($_,38,8);
		                	my $y_val = sprintf("%8.3f",$y);
		                	push @ycord1,0.1*$y_val;
		                	my $z=substr($_,46,8);
		                	my $z_val = sprintf("%8.3f",$z); 
		                	push @zcord1,0.1*$z_val;
		                	$i = $i+1;
		               
		                	my $res = substr($_,17,3);
		                	push @res_name,$res;
	 			
		         	}


	       		}
		}

		reset 'xcord1';
		reset 'ycord1';
		reset 'zcord1';
		close FH;
		$atom_type_file = ($ARGV[0]).'_atom_type.txt';
		open ( FG, "< $atom_type_file") or die "cannot open the file $!";
		my @at;
		while(<FG>)
		{
			my $a = substr($_,0,3);
			push @at,$a;

		}
		close FG;

		$len = scalar(@at);

		open (F1, "> $grofile")or die "cannot open the file $!";
		print F1 "$proteinname\n  $len\n";
		for($i = 0; $i<$len ; $i+=1)
		{
			printf F1 "%5d%5s%5s%5d%8.3f%8.3f%8.3f\n" , $i+1,$res_name[$i],$at[$i],$i+1,$xcord1[$i],$ycord1[$i],$zcord1[$i];
		}
		print F1 "  0.00000   0.00000   0.00000";

		close F1;

	}


	
	system("rm $trac_res $atom_type_file");
	system("rm ./exclude_res.txt");

	if($metal_flg == 1 && $metal_atom eq 'CA')
	{ system("rm $atom_type_file fort.2 m.txt metal.txt "); }

	if($metal_flg == 1 && $metal_atom eq 'CU')
	{ system("rm $atom_type_file m.txt metal.txt "); }

	if($metal_flg == 1 && $metal_atom eq 'ZN')
	{ system("rm $atom_type_file fort.2 m.txt metal.txt "); }

	if($metal_flg == 1 && $metal_atom eq 'FE')
	{ system("rm $atom_type_file  m.txt metal.txt "); }

	if($metal_flg == 1 && $metal_atom eq 'MG')
	{ system("rm $atom_type_file m.txt metal.txt "); }

	


}
print("\n## FINISH ##\n\n");
