#!/usr/bin/perl
use lib '.';
use elmo;
use Shell;
use File::Copy;
use File::Path;
use Getopt::Long;

our $VERSION = "2.1";

# VARIABLES TO EDIT
######## PUT THE GOOD PATH : Absolute path of the CHARMM program 
our $charmm = "/net/bccs/cbu/linux/Charmm/c33b1/exec/gnu/charmm";
#our $charmm = "/home/lbm/progmgr/charmm/c33b2.ia64/exec/gnu.scalar-intel/charmm";
#our $charmm = "/Users/annick/Prog/c32a1/exec/osx/charmm";
# Command line to execute CHARMM on mutiple processors (used by default for the execution of the molecular dynamics and the processing of the conformations)
 # To determine the number of processors and the nodes to be used, the template file tmp/mymachines.tmp must be edited.
our $charmm_mpi = "./mpirun -machinefile ./mymachines.LINUX -np 8 /usr/local/bin/charmm";
######## PUT THE GOOD PATH : Absolute path of the UHBD program
our $Uhbd = "/home/fimm/cbu/cgr061/pk_run_nov2005/pk/scripts";
#our $Uhbd = "/home/elyette/programmes/pk_run_nov2005/pk/scripts";
 # Array refering to the ligands of interest (not used for this version)
our @Hetatm = qw//;


# ANALYSIS OF THE OPTIONAL PARAMETERS
# See the documentation to obtain more detail for each option
my $Verbose = 1;
my $All = -1;
my $Base = -1;
my $Calc = -1;
my $Dyna = -1;
my $Conf = -1;
my $Int = 0;
my $Run = -1;
GetOptions (	'Debug' => sub { warn "\t\/!\\ DEBUG MODE \/!\\"; $| = 1 },
		'Verbose!' => \$Verbose,
		'Int!' => \$Int,
		'All!' => \$All,
			'Base!' => \$Base,
				'BaseBuild!' => \$BaseBuild,
				'BaseMiniHis!'  => \$BaseMiniHis,
				'BaseCrd!'   => \$BaseCrd,
			'SkipHis!' => \$SkipHis,
			'His!' => \$His,
                        'zapthat!' => \$zapthat,
			'Calc!' => \$Calc,
				'CalcMini!' => \$CalcMini,
				'CalcSas!'  => \$CalcSas,
				'CalcVdw!'  => \$CalcVdw,
				'CalcUhbd!' => \$CalcUhbd,
				'CalcGbmv!' => \$CalcGbmv,
				'CalcSumUhbd!'  => \$CalcSumUhbd,
				'CalcSumGbmv!'  => \$CalcSumGbmv,
			'FirstPlot' => \$FirstPlot,
			'Dyna!' => \$Dyna,
				'DynaMini!' => \$DynaMini,
				'DynaInp!'  => \$DynaInp,
				'DynaExec!' => \$DynaExec,
			'Conf!' => \$Conf,
				'ConfEner!' => \$ConfEner,
				'ConfTraj!' => \$ConfTraj,
				'Separated!' => \$Separated,
				'Together!' => \$Together,
				'ConfSas!' => \$ConfSas,
				'ConfVdw!' => \$ConfVdw,
				'ConfUHBD!' => \$ConfUhbd,
				'ConfGBMV!' => \$ConfGbmv,
			'Run!' => \$Run,
			'SumConf!' => \$SumConf,
			'average' => \$average,
			'stdev' => \$stdev,
			'supplement' => \$supplement,
			'average_nonpond' => \$average_nonpond,
			'ConfPlot' => \$ConfPlot,
			'mpi!' => \$mpi,
		'File=s' => \$filin,
		'Data=s' => \$data,
		'Complex=s' => \$Complex	);
if ($All == 1) {
	$BaseBuild = 1 unless (defined $BaseBuild); $BaseMiniHis = 1 unless (defined $BaseMiniHis); $BaseCrd = 1 unless (defined $BaseCrd);
	$His = 1 unless (defined $His);
	$CalcMini = 1 unless (defined $CalcMini); $CalcSas = 1 unless (defined $CalcSas); $CalcVdw = 1 unless (defined $CalcVdw); $CalcUhbd = 1 unless (defined $CalcUhbd); $CalcGbmv = 1 unless (defined $CalcGbmv); $CalcSumUhbd = 1 unless (defined $CalcSumUhbd); $CalcSumGbmv = 1 unless (defined $CalcSumGbmv);
	$ConfEner = 1 unless (defined $ConfEner); $ConfTraj = 1 unless (defined $ConfTraj); $ConfSas = 1 unless (defined $ConfSas); $ConfVdw = 1 unless (defined $ConfVdw);$ConfUhbd = 1 unless (defined $ConfUhbd); $ConfGbmv = 1 unless (defined $ConfGbmv);
	$DynaMini = 1 unless (defined $DynaMini); $DynaInp = 1 unless (defined $DynaInp); $DynaExec = 1 unless (defined $DynaExec);
	$SumConf = 1 unless (defined $SumConf); $Run = 1 unless (defined $Run);
}
elsif ($All == 0) {
	$BaseBuild = 0 unless (defined $BaseBuild); $BaseMini = 0 unless (defined $BaseMini); $BaseCrd = 0 unless (defined $BaseCrd);
	$His = 0 unless (defined $His);
	$CalcMini = 0 unless (defined $CalcMini); $CalcSas = 0 unless (defined $CalcSas); $CalcVdw = 0 unless (defined $CalcVdw); $CalcUhbd = 0 unless (defined $CalcUhbd); $CalcGbmv = 0 unless (defined $CalcGbmv); $CalcSumUhbd = 0 unless (defined $CalcSumUhbd); $CalcSumGbmv = 0 unless (defined $CalcSumGbmv);
	$ConfEner = 0 unless (defined $ConfEner); $ConfTraj = 0 unless (defined $ConfTraj);$ConfSas = 0 unless (defined $ConfSas); $ConfVdw = 0 unless (defined $ConfVdw); $ConfUhbd = 0 unless (defined $ConfUhbd); $ConfGbmv = 0 unless (defined $ConfGbmv);
	$DynaMini = 0 unless (defined $DynaMini); $DynaInp = 0 unless (defined $DynaInp); $DynaExec = 0 unless (defined $DynaExec);
	$SumConf = 0 unless (defined $SumConf); $Run = 0 unless (defined $Run);
}
if ($Base == 1) {
	$BaseBuild = 1 unless (defined $BaseBuild); $BaseMiniHis = 1 unless (defined $BaseMiniHis); $BaseCrd = 1 unless (defined $BaseCrd);
}
elsif ($Base == 0) {
	$BaseBuild = 0 unless (defined $BaseBuild); $BaseMiniHis = 0 unless (defined $BaseMiniHis); $BaseCrd = 0 unless (defined $BaseCrd);
}
if ($Calc == 1) {
	$CalcMini = 1 unless (defined $CalcMini); $CalcSas = 1 unless (defined $CalcSas); $CalcVdw = 1 unless (defined $CalcVdw); $CalcUhbd = 1 unless (defined $CalcUhbd); $CalcGbmv = 1 unless (defined $CalcGbmv); $CalcSumUhbd = 1 unless (defined $CalcSumUhbd); $CalcSumGbmv = 1 unless (defined $CalcSumGbmv);
}
elsif ($Calc == 0) {
	$CalcMini = 0 unless (defined $CalcMini); $CalcSas = 0 unless (defined $CalcSas); $CalcVdw = 0 unless (defined $CalcVdw); $CalcUhbd = 0 unless (defined $CalcUhbd); $CalcGbmv = 0 unless (defined $CalcGbmv); $CalcSumUhbd = 0 unless (defined $CalcSumUhbd); $CalcSumGbmv = 0 unless (defined $CalcSumGbmv);
}
if ($Conf == 1) {
	$ConfEner = 1 unless (defined $ConfEner); $ConfTraj = 1 unless (defined $ConfTraj); $ConfUhbd = 1 unless (defined
	$ConfUhbd); $ConfGBMV = 1 unless (defined $ConfGBMV);$ConfSas = 1 unless (defined $ConfSas); $ConfVdw = 1 unless (defined $ConfVdw);
}
elsif ($Conf == 0) {
	$ConfEner = 0 unless (defined $ConfEner); $ConfTraj = 0 unless (defined $ConfTraj); $ConfUhbd = 0 unless (defined
	$ConfUhbd); $ConfGbmv = 0 unless (defined $ConfGbmv);$ConfSas = 0 unless (defined $ConfSas); $ConfVdw = 0 unless (defined $ConfVdw);
}
if ($Dyna == 1) {
	$DynaMini = 1 unless (defined $DynaMini); $DynaInp = 1 unless (defined $DynaInp); $DynaExec = 1 unless (defined $DynaExec);
}
elsif ($Dyna == 0) {
	$DynaMini = 0 unless (defined $DynaMini); $DynaInp = 0 unless (defined $DynaInp); $DynaExec = 0 unless (defined $DynaExec);
}

if ($mpi==1) {
	$charmm =$charmm_mpi;
}


if ((!defined $filin)  && ($Int == 0)){
	print "Name of the input PDB file : "; chomp($filin = <STDIN>);
}
if (($CalcSumUhbd == 1) || ($CalcSumGbmv) || ($SumConf == 1)){
	print "Enter the name of the first residue of the second part of the system followed by a comma and the number of this residue :";
        $ref=$zapthat;#<STDIN>;
        @ref=split(',',$ref);
}
$start = time();
chomp (my $path = ".");
chomp (my $path2 = pwd());

# CREATION OF THE WORKING DIRECTORIES
if ($Int == 0) {
	mkpath (['inp','lib','out','scr'],0,0751);
	copy ("tmpsep/mini_sd.str","inp/mini_sd.str");
	copy ("tmpsep/mini_abnr.str","inp/mini_abnr.str");
	#copy ("tmpsep/allatom.str","inp/allatom.str");
	copy ("tmpsep/selec.def","lib/selec.def");
	#copy ("tmpsep/par_all22_prot_na.inp","lib/par_all22_prot_na.inp");
	#copy ("tmpsep/top_all22_prot_na.inp","lib/top_all22_prot_na.inp");
	#copy ("tmpsep/mymachines.tmp","mymachines.mpi");	
# OPENING OF THE INPUT PDB FILE
	open (FILIN, "< $filin") || die "Unable to open the file $filin : $!\n";
	chomp (my $path = pwd());
	my ($fname,$ext) = split ('\.',lc $filin);
# CREATION OF THE LOG FILE
	open (WARN, "> out/$fname.out") || warn "Unable to create the file out/$fname.out.\n";
	printf (WARN " --------------------------------------------------------------------\n");
	printf (WARN "| Log from the execution of the script Leto.pl on the file %s |\n",$filin);
	printf (WARN " --------------------------------------------------------------------\n");
	warn "\nProcessing the file $filin in $path ...\n" if ($Verbose == 1);	
	my @lines = <FILIN>;
	close (FILIN);

# READING OF THE FILE
	warn "\nReading of the file :\n" if ($Verbose == 1);
	$i=0;
	while ($lines[$i]!~/^ATOM/) {
		$i++;
	}
	if ((substr($lines[$i],21,1)=~/[a-zA-Z0-9]/)) { # for the case of pdb from the PDB
		print "\t\tPDB\n";
		($ref_chains, $ref_res, $ref_atoms, $ref_last, $hetatm, $nb_ch_sys, $nb_het_sys, $water) = Leto::GrpAtoms->read_pdb(\@lines,$Complex);
		$type=1; # type one for the PDB
	}
	else  { # for the case of pdb from charmm or modeller
		print "\t\tcharmm/modeller\n";
		($ref_chains, $ref_res, $ref_atoms, $ref_last, $hetatm, $nb_ch_sys,$nb_het_sys, $water) = Leto::GrpAtoms->read_pdb_charmm(\@lines,$Complex);
		$type=2; # type two for the charmm_pdb or modeller_pdb
	}
	printf (WARN "\n>> File successfully read.\n");
	foreach $line (0..$#{$hetatm}) {
		$hetatm->[$line] =~ s/\s+/ /g;
		printf (WARN "\tHetero-Atom deleted : %s\n", $hetatm->[$line]);
	}

# CORRECTION OF THE FILE TO COMPLY WITH THE CHARMM STANDARD
	warn "\nCorrection of file to comply with the CHARMM standard ...\n" if ($Verbose == 1);
	($ref_atoms,$before,$after) = Leto::GrpAtoms->pdb2charmm($ref_atoms,$ref_last,$type);
	printf (WARN "\n>> File corrected to be used in CHARMM : \n\tHIS -> HSP\n\tCD1 ILE -> CD ILE\n\tOXT -> OT2\n\tO TER -> OT1\n");
	printf (WARN "> If an atom has been recorded in the PDB file with more than 1 set of coordinates,\n\t the program determines and selects the set of coordinates with the highest occupancy.\n\tBy default for the same occupancy of a given atom,\n\t the first set encountered in the PDB file is selected.\n\tThe Atoms objects from the other sets of coordinates are just deleted.\n");
	foreach $step (0..$#{$after}) {
		if (($before->[$step] ne $after->[$step]) && ($before->[$step] ne $before->[$step+1])) {
			$before->[$step] =~ s/\s+/ /g; $after->[$step] =~ s/\s+/ /g;
			if ($after->[$step] =~ /HSP/) {
				printf (WARN "\tResidue  %s modified in  %s\n", $before->[$step], $after->[$step]);
			}
			elsif ($after->[$step] =~ /alt_location/) {
				printf (WARN "\tAtom %s had alternative locations\n", $before->[$step]);
			}
			else {
				printf (WARN "\tAtom %s modified in %s\n", $before->[$step], $after->[$step]);
			}
		}
	}

# RENUMBERING OF THE FILE FROM 1
	warn "\nRenumbering of the file ...\n" if ($Verbose == 1);
	($ref_atoms,$ref_last) = Leto::GrpAtoms->renum_pdb($ref_atoms,$ref_last,$Verbose,$type);
	printf (WARN "\n>> File with atoms and residues renumbered from 1.\n");

# SPLIT AND WRITING OF SEVERAL OUTPUT PDB FILES
	warn "\nAutomatic generation of the output PDB files ...\n" if ($Verbose == 1);
	$ref_pdb = Leto::GrpAtoms->split_pdb2charmm($ref_atoms,$filin,$type);
	printf (WARN "\n>> Output PDB files automatically generated.\n");
	if ($Verbose == 1) {
		printf ("\t%d PDB files were created in the scr/ directory : \n", $#{$ref_pdb}+1);
	}
	printf (WARN "\t%d PDB files were created in the scr/ directory : ", $#{$ref_pdb}+1);
	foreach $p (0..$#{$ref_pdb}) {
		printf (WARN ", %s.pdb",$ref_pdb->[$p]);
		printf ("\t\t------> %s.pdb\n",$ref_pdb->[$p]);
		}
	printf (WARN ".\n");
	printf (WARN "\tThe file %s.pdb contains the coordinates of the water molecules renamed after TIP3 residue name to comply with the CHARMM standard.\n", $ref_pdb->[-1]);

# Creation of the complex PDB file name : 
	foreach $i (0..$nb_ch_sys-1) {
		($pdb,$concatenation) = split('\_',$ref_pdb->[$i]);
		$com=$com . $concatenation;
		$ref=$i;
	}
	if ($nb_het_sys !=0) {
		$ref++;
		($pdb,$het) = split('\_',$ref_pdb->[$ref]);
		$com=$com . $het;
	}
	($pdb,$concatenation) = split('\_',$ref_pdb->[0]);
	push (@{$ref_pdb}, join('_',$pdb,$com));
	($pdb_name,)=split(/_/,$ref_pdb->[1]);
		
	print "\nEnter the name of the topology file (it must be in a lib directory) [default top_all22_prot.inp]: ";
	chomp ($top = <STDIN>);
	if ($top eq "") {
		$top="top_all22_prot.inp";
	}
	print "Enter the name of the parameter file (it must be in a lib directory) [default par_all22_prot.inp]: ";
	chomp ($par = <STDIN>);
	if ($par eq "") {
		$par="par_all22_prot.inp";
	}
	if (($BaseCrd ==1) || ($CalcMini == 1) || ($CalcSas == 1) || ($CalcVdw == 1) || ($CalcGbmv == 1)) {
		# asks to the user to define with charmm keyword each part of his system
		print "Enter the definition with charmm keyword of the first part :";
		chomp ($fp=<STDIN>);
		print "Enter the definition with charmm keyword of the second part :";
		chomp ($sp=<STDIN>);
	}	

# Reading of the ter file if it exists
	if (-e "dfsep/ter.dat") {
		open (TER, "< dfsep/ter.dat") || die "Unable to open the file dfsep/ter.dat : $!\n";	
		chomp(@ter=<TER>);
		close(TER);	
	}
	# Creating of a table containing the chain, first patch and last patch
	@term;
	foreach $i(1..$#ter) {
		($term[$i][0],$term[$i][1],$term[$i][2])=split(/\t/,$ter[$i]);
		$term[$i][0]=uc($term[$i][0]);
		$term[$i][1]=uc($term[$i][1]);
		$term[$i][2]=uc($term[$i][2]);
	}
	
	# If the patch for a chain is not in the ter.dat file, default patchs are NTER and CTER
	foreach $i (0..$nb_ch_sys-1) {
		($pdb,$concatenation) = split('\_',$ref_pdb->[$i]);
		foreach $a(1..$#term) {
			$ok=0;
			if (uc($term[$a][0]) eq uc($concatenation)) {
				$ok=1;
				last;
			}
		}
		if ($ok == 0) {
			$index=$#term+1;
			$term[$index][0]=uc($concatenation);
			$term[$index][1]=""; # used for print_build
			$term[$index][2]="";
			$term[$index][3]="NTER"; # used in his file for crd2uhbd
                        $term[$index][4]="CTER";

		}
	}
	$term_ref=\@term;
	
# Reading of the patchs file if it exists
	if (-e "dfsep/patch.dat") {
		open (PATCH, "< dfsep/patch.dat") || die "Unable to open the file dfsep/patch.dat : $!\n";	
		chomp(@patchs=<PATCH>);
		close(PATCH);	
	}
	# Creating of a table containing the chain, name of the patch and residu number
	@patch;
	foreach $j(1..$#patchs) {
		($patch[$j][0],$patch[$j][1],$patch[$j][2])=split(/\t/,$patchs[$j]);
	}
	$patch_ref=\@patch;
}

# CREATION OF THE INPUT BUILD FILE FROM TEMPLATE
if ($BaseBuild == 1) {
	warn "\nGeneration of the file build.inp from the template to proceed to the Hbuild module of CHARMM ...\n" if ($Verbose == 1);
	$build="base"; # this variable allows to put the good name in the print_build ($chain.pdb) to not confuse with the new build after that protonation state of his were know
	($ref_atoms,$ref_pdb) = Leto::GrpAtoms->print_build($ref_atoms,$ref_pdb,$path,$Verbose,$nb_ch_sys,$nb_het_sys,$water,$type,$build,$pdb_name,$top,$par,$term_ref,$patch_ref);
	printf (WARN "\n>> Input Hbuild CHARMM module inp/$pdb_name\_complex.build.inp created from template.\n");
# EXECUTION OF THE HBUILD CHARMM MODULE
	$begin = time();
	if ($Verbose == 1) {
		warn "\nNow building the hydrogen atoms using the Hbuild CHARMM module ...\n";
		warn "\tcharmm < inp/$pdb_name\_complex.build.inp > out/$pdb_name\_complex.build.out\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
	}
	system ("$charmm < inp/$pdb_name\_complex.build.inp >out/$pdb_name\_complex.build.out");
	$end = time();
	$sec = ($end-$begin) %60;
	$dif = (($end-$begin)-$sec) /60;
	$min = $dif %60;
	($result,$ref_lines) = Leto::GrpAtoms->out_verif("$pdb_name\_complex.build",$path,$Verbose);
	foreach $l (0..$#{$ref_lines}) {
		printf (WARN "%s", $ref_lines->[$l]);
	}
	if ($result == 1) {
		warn "$ref_lines->[0]" if ($Verbose == 1);
	}
	else {
		die "\t/!\\$ref_lines->[0]";
	}
	$ref_lines = undef;
	printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min, $sec) if ($Verbose == 1);
	warn "\tSee the log : out/$pdb_name\_complex.build.out.\n\n" if ($Verbose == 1);
	printf(WARN "\n>> Execution of the Hbuild CHARMM module completed in %d minutes and %d seconds.\n\tSee the log : out/%s.build.out.\n\n", $min, $sec, $pdb_name);
}


# CREATION OF THE INPUT SMALL MINIMIZATION FILE BEFORE PROTONATION STATE OF HIS CALCULATION
if ($BaseMiniHis == 1) {
	warn "\nGeneration of the file mini.inp from the template to proceed to the minimization with CHARMM ...\n" if ($Verbose == 1);
	$mini_ext = "small";
	Leto::GrpAtoms->print_mini($mini_ext,$path,$Verbose,$nb_ch_sys,$nb_het_sys,$pdb_name,$top,$par);
	printf (WARN "\n>> Input minimization file inp/$pdb_name\_complex.mini_$mini_ext.inp created from template.\n") if ($Verbose == 1);
# EXECUTION OF THE MINIMIZATION WITH CHARMM
	$begin = time();
	if ($Verbose == 1) {
		warn "\nNow minimizing with CHARMM ...\n";
		warn "\tcharmm < inp/$pdb_name\_complex.mini_$mini_ext.inp > out/$pdb_name\_complex.mini_$mini_ext.out\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
	}
	system ("$charmm < inp/$pdb_name\_complex.mini_$mini_ext.inp > out/$pdb_name\_complex.mini_$mini_ext.out");
	$end = time();
	$sec = ($end-$begin) %60;
	$dif = (($end-$begin)-$sec) /60;
	$min = $dif %60;
	($result,$ref_lines) = Leto::GrpAtoms->out_verif("$pdb_name\_complex.mini_$mini_ext",$path,$Verbose);
	foreach $l (0..$#{$ref_lines}) {
		printf (WARN "%s", $ref_lines->[$l]);
	}
	if ($result == 1) {
		warn "$ref_lines->[0]" if ($Verbose == 1);
	}
	else {
		die "\t/!\\$ref_lines->[0]";
	}
	$ref_lines = undef;
	printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min, $sec) if ($Verbose == 1);
	warn "\tSee the log : out/$pdb_name\_complex.mini_$mini_ext.out.\n\n" if ($Verbose == 1);
	printf (WARN "\n>> Minimization with CHARMM completed in %d minutes and %d seconds.\n\tSee the log : out/%s.mini_$mini_ext.out.\n\n", $min, $sec, $pdb_name) if ($Verbose == 1);
}


# CREATION OF THE INPUT CHARGES-RADII WRITING FILE
if ($BaseCrd == 1) {
	# Check if a $pdb_name.mini_small.crd exists, that means the option baseminihis was used
	# Allows to know if this file exist when the option baseminihis and basecrd are used separatly
	if (-e "scr/$pdb_name\_complex.mini_small.crd") {
		$mini_ext="small";
	}
	warn "\nGeneration of the file write_cha_rad.inp from template to proceed to the charges-radii writing with CHARMM ...\n" if ($Verbose == 1);
	Leto::GrpAtoms->print_crd($mini_ext,$path,$Verbose,$nb_ch_sys,$nb_het_sys,$pdb_name,$top,$par,$fp,$sp);
	printf (WARN "\n>> Input charges-radii writing file inp/$pdb_name\_complex.crd.inp created from template.\n");
# PROCESSING THE CHARGES-RADII WRITING WITH CHARMM
	$begin = time();
	if ($Verbose == 1) {
		warn "\nNow writing the charges-radii CRD files with CHARMM ...\n";
		warn "\tcharmm < inp/$pdb_name\_complex.crd.inp > out/$pdb_name\_complex.crd.out\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
	}
	system ("$charmm < inp/$pdb_name\_complex.crd.inp > out/$pdb_name\_complex.crd.out");
	$end = time();
	$sec = ($end-$begin) %60;
	$dif = (($end-$begin)-$sec) /60;
	$min = $dif %60;
	($result,$ref_lines) = Leto::GrpAtoms->out_verif("$pdb_name\_complex.crd",$path,$Verbose);
	foreach $l (0..$#{$ref_lines}) {
		printf (WARN "%s", $ref_lines->[$l]);
	}
	if ($result == 1) {
		warn "$ref_lines->[0]" if ($Verbose == 1);
	}
	else {
		die "\t/!\\$ref_lines->[0]";
	}
	$ref_lines = undef;
	printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min,$sec) if ($Verbose == 1);
	warn "\tSee the log : out/$pdb_name\_complex.crd.out.\n\n" if ($Verbose == 1);
	printf (WARN "\n>> Writing of the CRD files with CHARMM completed in %d seconds.\n\tSee the log : out/%s.crd.out.\n\n", $sec, $pdb_name);
}


# PATCHING THE CRD FILES TO COMPLY WITH THE UHBD STANDARD WITHOUT MAKE UHBD CALCULATION FOR THE HISTIDINE
# Use this option when you know already the protonation state of histidine  
if ($SkipHis  == 1) {
	if (-e "scr/$pdb_name\_complex.mini_small.crd") {
		$mini_ext="small";
	}
	warn "\nCorrection of the CRD files to comply with the UHBD standard ...\n" if ($Verbose == 1);
	($before,$after) = Leto::GrpAtoms->crd2uhbd($ref_atoms,$ref_last,"scr/*-{cha,rad}.crd",$term_ref,$patch_ref);
	printf (WARN "\n>> CRD files corrected to be used in UHBD.\n");
	printf (WARN "> Backbone atoms of the first and the last residues of each chain renamed in NTER or CTER.\n");
	printf (WARN "> Atoms of the cysteins involved in disulfide bonds renamed in CSS.\n");
	foreach $step (0..$#{$after}) {
		if (($before->[$step] ne $after->[$step]) && ($before->[$step] ne $before->[$step+1])) {
			$before->[$step] =~ s/\s+/ /g; $after->[$step] =~ s/\s+/ /g;
			if ($after->[$step] =~ /CSS/) {
				printf (WARN "\tResidue %s modified in %s\n", $before->[$step], $after->[$step]);
			}
			else {	
				printf (WARN "\tAtom  %s modified in  %s\n", $before->[$step], $after->[$step]);
			}
		}
	}
	# READING OF THE RESULT FILES AND MODIFICATION OF THE HISTIDINES IN CONSEQUENCE
	warn "\nCorrection of each histidine according to its protonation state ...\n" if ($Verbose == 1);
	($his,$before,$after) = Leto::GrpAtoms->his_skip($ref_pdb,$pdb_name,$mini_ext);
	printf (WARN "\n>> Histidines corrected according to their protonation state.\n");
	printf (WARN "> The superfluous hydrogen atoms were deleted.\n");
	foreach $step (0..$#{$after}) {
		if ($before->[$step] =~ /Chain /) {
		printf (WARN "\t%s has no Histidine.\n\tResult files skipped.\n", uc($after->[$step]));
		}
		elsif ($before->[$step] =~ /(Pk|Xave)/) {			
			my ($file,$chain) = split(' ',$after->[$step]);
			printf (WARN "\tNo %s result file found for chain %s.\n\tStep skipped.\n", $file, uc($chain));
		}
		else {
			$before->[$step] =~ s/\s+/ /g; $after->[$step] =~ s/\s+/ /g;
			printf (WARN "\tResidue  %s modified\n", $after->[$step]);
			printf (WARN "\tAtom  %s deleted from  %s\n", $before->[$step], $after->[$step]);
		}
	}
	for ($step=0; $step<=$#{$his}; $step++) {
		my ($pdb,$chain) = split("\_", $ref_pdb->[$his->[$step+1]]);
		printf (WARN "\tResidue %d has different protonation states : \n\t%s in chain %s and %s in the complex.\n", $his->[$step], $his->[$step+2], uc($chain), $his->[$step+3]);
		$step += 3;
	}
 	
	#CREATION OF THE CRD AND PSF FILES WITH THE HISTIDINES IN THE CORRECT PROTONATION STATE
	$begin = time();
	warn "\nGeneration of the file build_psf_crd.inp from the template to create the CRD and PSF files with CHARMM ...\n" if ($Verbose == 1);
	$build="his"; # this variable allows to put the good name in the print_build ($chain-cha.crd.his) to not confuse with the first build
	($ref_atoms,$ref_pdb) = Leto::GrpAtoms->print_build($ref_atoms,$ref_pdb,$path,$Verbose,$nb_ch_sys,$nb_het_sys,$water,$type,$build,$pdb_name,$top,$par,$term_ref,$patch_ref);
	printf (WARN "\n>> Input build file inp/$pdb_name\_complex.build_psf_crd.inp created from template.\n") if ($Verbose == 1);
	if ($Verbose == 1) {
		warn "\nNow building CRD and PSF files with its histidines in their correct protonation state with CHARMM ...\n";
		warn "\tcharmm < inp/$pdb_name\_complex.build_psf_crd.inp > out/$pdb_name\_complex.build_psf_crd.out\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
	}
	system ("$charmm < inp/$pdb_name\_complex.build_psf_crd.inp > out/$pdb_name\_complex.build_psf_crd.out");
	$end = time();
	$sec = ($end-$begin) %60;
	$dif = (($end-$begin)-$sec) /60;
	$min = $dif %60;
	($result,$ref_lines) = Leto::GrpAtoms->out_verif("$pdb_name\_complex.build_psf_crd",$path,$Verbose);
	foreach $l (0..$#{$ref_lines}) {
		printf (WARN "%s", $ref_lines->[$l]);
	}
	if ($result == 1) {
		warn "$ref_lines->[0]" if ($Verbose == 1);
	}
	else {
		die "\t/!\\$ref_lines->[0]";
	}
	$ref_lines = undef;
	printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min,$sec) if ($Verbose == 1);
	warn "\tSee the log : out/$pdb_name\_complex.build_psf_crd.out.\n\n" if ($Verbose == 1);
	printf (WARN "\n>> Writing of the CRD and PSF files with the histidines in the correct protonation state with CHARMM completed in %d seconds.\n\tSee the log : out/%s.build_psf_crd.out.\n\n", $sec, $ref_pdb->[-1]);
}


# PATCHING THE CRD FILES TO COMPLY WITH THE UHBD STANDARD
if ($His == 1) {
	 # Check if a $pdb_name.mini_small.crd exists, that means the option baseminihis was used
        # Allows to know if this file exist when the option baseminihis and basecrd are used separatly
        if (-e "scr/$pdb_name\_complex.mini_small.crd") {
                $mini_ext="small";
        }
	warn "\nCorrection of the CRD files to comply with the UHBD standard ...\n" if ($Verbose == 1);
	($before,$after) = Leto::GrpAtoms->crd2uhbd($ref_atoms,$ref_last,"scr/*-{cha,rad}.crd",$term_ref,$pacth_ref);
	printf (WARN "\n>> CRD files corrected to be used in UHBD.\n");
	printf (WARN "> Backbone atoms of the first and the last residues of each chain renamed in NTER or CTER.\n");
	printf (WARN "> Atoms of the cysteins involved in disulfide bonds renamed in CSS.\n");
	foreach $step (0..$#{$after}) {
		if (($before->[$step] ne $after->[$step]) && ($before->[$step] ne $before->[$step+1])) {
			$before->[$step] =~ s/\s+/ /g; $after->[$step] =~ s/\s+/ /g;
			if ($after->[$step] =~ /CSS/) {
				printf (WARN "\tResidue %s modified in %s\n", $before->[$step], $after->[$step]);
			}
			else {
				printf (WARN "\tAtom  %s modified in  %s\n", $before->[$step], $after->[$step]);
			}
		}
	}
	
	# CREATION OF THE RUN.COM FILE TO PROCEED TO UHBD
	my $boucle;
	$grid = 0.3; #the initial grid value.
	while ($boucle == 0) { #loop allowing to remake script run.com and calculation if the initial grid were too small.
	
		warn "\nGeneration of the file run.com from template to process with UHBD on the different chains with grid=$grid...\n" if ($Verbose == 1);
		Leto::GrpAtoms->print_run($Uhbd,$Verbose,$grid,$pdb_name);
		printf (WARN "\n>> UHBD running command file inp/$pdb_name\_complex.run.com created from template.\n");
	
	# PROCESSING UHBD
		$begin = time();
		foreach $p ("fp","sp","cx") {
			
			mkpath ("uhbd/his/$pdb_name\_$p/tmp",0,0751);
			if ($mini_ext eq "small") {
				copy ("scr/$pdb_name-mini\_$p-cha.crd","uhbd/his/$pdb_name\_$p/$pdb_name\_$p-cha.crd");
				copy ("scr/$pdb_name-mini\_$p-rad.crd","uhbd/his/$pdb_name\_$p/$pdb_name\_$p-rad.crd");
			}
			else {
				copy ("scr/$pdb_name\_$p-cha.crd","uhbd/his/$pdb_name\_$p/$pdb_name\_$p-cha.crd");
				copy ("scr/$pdb_name\_$p-rad.crd","uhbd/his/$pdb_name\_$p/$pdb_name\_$p-rad.crd");
			}
			

		}
		if ($Verbose == 1) {
		warn "\nNow running UHBD to determine the protonation state of the Histidines ...\n";
		warn "\tcsh inp/$pdb_name\_complex.run.com > out/$pdb_name\_complex.run.out\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
		}
		system ("csh inp/$pdb_name\_complex.run.com > out/$pdb_name\_complex.run.out");
		$end = time();
		$sec = ($end-$begin) %60;
		$dif = (($end-$begin)-$sec) /60;
		$min = $dif %60;
		$dif = ($dif-$min) /60;
		$h = $dif %24;
		($result,$ref_lines) = Leto::GrpAtoms->out_verif("$pdb_name\_complex.run",$path,$Verbose);
		foreach $l (0..$#{$ref_lines}) {
			printf (WARN "%s", $ref_lines->[$l]);
		}
		if ($result == 1) { # calculation was finished without problem.
			warn "$ref_lines->[0]" if ($Verbose == 1);
			$boucle=1;
		}
		elsif ($result == 2) { # the grid is too small, the run.com is recreated with a value of grid adapted.
			$grid=substr($ref_lines->[1],43,5);
			$grid="0".$grid;
			foreach $p ("fp","sp","cx") {
				system "rm -r uhbd/his/$pdb_name\_$p/tmp/*";
			}
			$boucle=0;
		}
		else { # there was a problem during calculation.
			die "\t/!\\$ref_lines->[0]";
		}
	}
	printf ("Execution of UHBD completed in %d hours %d minutes and %d seconds.\n", $h, $min, $sec) if ($Verbose == 1);
	warn "\tSee the log : out/$pdb_name\_complex.run.out.\n\n" if ($Verbose == 1);
	printf (WARN "\n>> Execution of UHBD completed in %d hours %d minutes and %d seconds.\n\tSee the log : out/%s\_complex.run.out.\n\n", $h, $min, $sec,$pdb_name);

	# READING OF THE RESULT FILES AND MODIFICATION OF THE HISTIDINES IN CONSEQUENCE
	warn "\nCorrection of each histidine according to its protonation state ...\n" if ($Verbose == 1);
	($his,$before,$after) = Leto::GrpAtoms->his_uhbd($ref_pdb,$pdb_name,$mini_ext);
	printf (WARN "\n>> Histidines corrected according to their protonation state.\n");
	printf (WARN "> The superfluous hydrogen atoms were deleted.\n");
	foreach $step (0..$#{$after}) {
		if ($before->[$step] =~ /Chain /) {
		printf (WARN "\t%s has no Histidine.\n\tResult files skipped.\n", uc($after->[$step]));
		}
		elsif ($before->[$step] =~ /(Pk|Xave)/) {			
			my ($file,$chain) = split(' ',$after->[$step]);
			printf (WARN "\tNo %s result file found for chain %s.\n\tStep skipped.\n", $file, uc($chain));
		}
		else {
			$before->[$step] =~ s/\s+/ /g; $after->[$step] =~ s/\s+/ /g;
			printf (WARN "\tResidue  %s modified\n", $after->[$step]);
			printf (WARN "\tAtom  %s deleted from  %s\n", $before->[$step], $after->[$step]);
		}
	}
	for ($step=0; $step<=$#{$his}; $step++) {
		my ($pdb,$chain) = split("\_", $ref_pdb->[$his->[$step+1]]);
		printf (WARN "\tResidue %d has different protonation states : \n\t%s in chain %s and %s in the complex.\n", $his->[$step], $his->[$step+2], uc($chain), $his->[$step+3]);
		$step += 3;
	}
 	
	#CREATION OF THE CRD AND PSF FILES WITH THE HISTIDINES IN THE CORRECT PROTONATION STATE
	$begin = time();
	warn "\nGeneration of the file build_psf_crd.inp from the template to create the CRD and PSF files with CHARMM ...\n" if ($Verbose == 1);
	$build="his";
	($ref_atoms,$ref_pdb) = Leto::GrpAtoms->print_build($ref_atoms,$ref_pdb,$path,$Verbose,$nb_ch_sys,$nb_het_sys,$water,$type,$build,$pdb_name,$top,$par,$term_ref,$patch_ref);
	printf (WARN "\n>> Input build file inp/$pdb_name\_complex.build_psf_crd.inp created from template.\n") if ($Verbose == 1);
	if ($Verbose == 1) {
		warn "\nNow building CRD and PSF files with its histidines in their correct protonation state with CHARMM ...\n";
		warn "\tcharmm < inp/$pdb_name\_complex.build_psf_crd.inp > out/$pdb_name\_complex.build_psf_crd.out\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
	}
	system ("$charmm < inp/$pdb_name\_complex.build_psf_crd.inp > out/$pdb_name\_complex.build_psf_crd.out");
	$end = time();
	$sec = ($end-$begin) %60;
	$dif = (($end-$begin)-$sec) /60;
	$min = $dif %60;
	($result,$ref_lines) = Leto::GrpAtoms->out_verif("$pdb_name\_complex.build_psf_crd",$path,$Verbose);
	foreach $l (0..$#{$ref_lines}) {
		printf (WARN "%s", $ref_lines->[$l]);
	}
	if ($result == 1) {
		warn "$ref_lines->[0]" if ($Verbose == 1);
	}
	else {
		die "\t/!\\$ref_lines->[0]";
	}
	$ref_lines = undef;
	printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min,$sec) if ($Verbose == 1);
	warn "\tSee the log : out/$pdb_name\_complex.build_psf_crd.out.\n\n" if ($Verbose == 1);
	printf (WARN "\n>> Writing of the CRD and PSF files with the histidines in the correct protonation state with CHARMM completed in %d seconds.\n\tSee the log : out/%s.build_psf_crd.out.\n\n", $sec, $ref_pdb->[-1]);
}


# CREATION OF THE INPUT MINIMIZATION FILE BEFORE THE ENERGY CALCULATION
if ($CalcMini == 1) {
	mkpath (['scr/mini'],0,0751);
	copy ("tmpsep/mini_sd.str","inp/mini_sd.str");
	copy ("tmpsep/mini_abnr.str","inp/mini_abnr.str");
	# CREATION OF INPUT FOR MINIMIZATION
	warn "\nGeneration of the file mini.inp from the template to proceed to the minimization with CHARMM ...\n" if ($Verbose == 1);
	$mini_ext = "his";
	Leto::GrpAtoms->print_mini($mini_ext,$path,$Verbose,$nb_ch_sys,$nb_het_sys,$pdb_name,$top,$par);
	printf (WARN "\n>> Input minimization file inp/$pdb_name\_complex.mini_$mini_ext.inp created from template.\n") if ($Verbose == 1);
	# EXECUTION OF THE MINIMIZATION WITH CHARMM
	$begin = time();
	if ($Verbose == 1) {
		warn "\nNow minimizing with CHARMM ...\n";
		warn "\tcharmm < inp/$pdb_name\_complex.mini_$mini_ext.inp > out/$pdb_name\_complex.mini_$mini_ext.out\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
	}
	system ("$charmm < inp/$pdb_name\_complex.mini_$mini_ext.inp > out/$pdb_name\_complex.mini_$mini_ext.out");
	$end = time();
	$sec = ($end-$begin) %60;
	$dif = (($end-$begin)-$sec) /60;
	$min = $dif %60;
	($result,$ref_lines) = Leto::GrpAtoms->out_verif("$pdb_name\_complex.mini_$mini_ext",$path,$Verbose);
	foreach $l (0..$#{$ref_lines}) {
		printf (WARN "%s", $ref_lines->[$l]);
	}
	if ($result == 1) {
		warn "$ref_lines->[0]" if ($Verbose == 1);
	}
	else {
		die "\t/!\\$ref_lines->[0]";
	}
	$ref_lines = undef;
	printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min, $sec) if ($Verbose == 1);
	warn "\tSee the log : out/$pdb_name\_complex.mini_$mini_ext.out.\n\n" if ($Verbose == 1);
	printf (WARN "\n>> Minimization with CHARMM completed in %d minutes and %d seconds.\n\tSee the log : out/%s.mini_$mini_ext.out.\n\n", $min, $sec, $ref_pdb->[-1]) if ($Verbose == 1);
	
	copy ("scr/mini/$pdb_name\_complex.mini_his_abnr_1000.crd","scr/$pdb_name\_complex.mini_his.crd");
	 
	# CREATION OF THE INPUT WHICH ALLOW TO CREATE CHA AND RAD FILES
	warn "\nGeneration of the file write_cha_rad.inp from template to proceed to the charges-radii writing with CHARMM ...\n" if ($Verbose == 1);
	Leto::GrpAtoms->print_crd($mini_ext,$path,$Verbose,$nb_ch_sys,$nb_het_sys,$pdb_name,$top,$par,$fp,$sp);
	printf (WARN "\n>> Input charges-radii writing file inp/$pdb_name\_complex.crd-mini.inp created from template.\n");
	# PROCESSING THE CHARGES-RADII WRITING WITH CHARMM
	$begin = time();
	if ($Verbose == 1) {
		warn "\nNow writing the charges-radii CRD files with CHARMM ...\n";
		warn "\tcharmm < inp/$pdb_name\_complex.crd-mini.inp > out/$pdb_name\_complex.crd-mini.out\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
	}
	system ("$charmm < inp/$pdb_name\_complex.crd-mini.inp > out/$pdb_name\_complex.crd-mini.out");
	$end = time();
	$sec = ($end-$begin) %60;
	$dif = (($end-$begin)-$sec) /60;
	$min = $dif %60;
	($result,$ref_lines) = Leto::GrpAtoms->out_verif("$pdb_name\_complex.crd-mini",$path,$Verbose);
	foreach $l (0..$#{$ref_lines}) {
		printf (WARN "%s", $ref_lines->[$l]);
	}
	if ($result == 1) {
		warn "$ref_lines->[0]" if ($Verbose == 1);
	}
	else {
		die "\t/!\\$ref_lines->[0]";
	}
	$ref_lines = undef;
	printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min,$sec) if ($Verbose == 1);
	warn "\tSee the log : out/$pdb_name\_complex.crd-mini.out.\n\n" if ($Verbose == 1);
	printf (WARN "\n>> Writing of the CRD files with CHARMM completed in %d seconds.\n\tSee the log : out/%s.crd-mini.out.\n\n", $sec, $pdb_name);
	
	# CORRECTION OF THE CRD
	warn "\nCorrection of the CRD files to comply with the UHBD standard ...\n" if ($Verbose == 1);
	($before,$after) = Leto::GrpAtoms->crd2uhbd($ref_atoms,$ref_last,"scr/*-mini_*-{cha,rad}.crd",$term_ref,$patch_ref);
	printf (WARN "\n>> CRD files corrected to be used in UHBD.\n");
	printf (WARN "> Backbone atoms of the first and the last residues of each chain renamed in NTER or CTER.\n");
	printf (WARN "> Atoms of the cysteins involved in disulfide bonds renamed in CSS.\n");
	foreach $step (0..$#{$after}) {
		if (($before->[$step] ne $after->[$step]) && ($before->[$step] ne $before->[$step+1])) {
			$before->[$step] =~ s/\s+/ /g; $after->[$step] =~ s/\s+/ /g;
			if ($after->[$step] =~ /CSS/) {
				printf (WARN "\tResidue %s modified in %s\n", $before->[$step], $after->[$step]);
			}
			else {
				printf (WARN "\tAtom  %s modified in  %s\n", $before->[$step], $after->[$step]);
			}
		}
	}
	
	# Allows to put the good names for the next steps
	foreach $p("fp","sp","cx") {
		rename ("scr/$pdb_name-mini\_$p-cha.crd","scr/$pdb_name-mini\_$p-cha.crd.his.ter");
		rename ("scr/$pdb_name-mini\_$p-rad.crd","scr/$pdb_name-mini\_$p-rad.crd.his.ter");
	}	
}


# CREATION OF THE INPUT SAS FILE
if ($CalcSas == 1) {
	if (-d "scr/mini") {
		$ext="mini";
	}
	warn "\nGeneration of the file sas.inp from template to process with the sas calculation with CHARMM ...\n" if ($Verbose == 1);
	Leto::GrpAtoms->print_sas($pdb_name,$path,$Verbose,$Int,$top,$par,$fp,$sp,$ext,"");
	printf (WARN "\n>> Input sas file inp/$pdb_name\_complex.sas.inp created from template.\n");
	
	# PROCESSING THE SAS CALCULATION WITH CHARMM
	$begin = time();
	if ($Verbose == 1) {
		warn "\nNow calculating the solvent accessible surface energy with CHARMM...\n";
		warn "\tcharmm < inp/$pdb_name\_complex.sas.inp > out/$pdb_name\_complex.sas.out\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
	}
	system ("$charmm < inp/$pdb_name\_complex.sas.inp > out/$pdb_name\_complex.sas.out");
	$end = time();
	$sec = ($end-$begin) %60;
	$dif = (($end-$begin)-$sec) /60;
	$min = $dif %60;
	($result,$ref_lines) = Leto::GrpAtoms->out_verif("$pdb_name\_complex.sas",$path,$Verbose);
	foreach $l (0..$#{$ref_lines}) {
		printf (WARN "%s", $ref_lines->[$l]);
	}
	if ($result == 1) {
		warn "$ref_lines->[0]" if ($Verbose == 1);
	}
	else {
		die "\t/!\\$ref_lines->[0]";
	}
	$ref_lines = undef;
	printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min, $sec) if ($Verbose == 1);
	warn "\tSee the log : out/$pdb_name\_complex.sas.out.\n\n" if ($Verbose == 1);
	printf (WARN "\n>> Calculation of the solvent accessible surface energy with CHARMM completed in %d minutes and %d seconds.\n\tSee the log : out/%s.sas.out.\n\n", $min, $sec, $ref_pdb->[-1]);
}


# CREATION OF THE INPUT VDW FILE
if ($CalcVdw == 1) {
	if (-d "scr/mini") {
		$ext="mini";
	}
	warn "\nGeneration of the file vdw.inp from template to process with the Van der Waals and electrostatic energies calculation with CHARMM ...\n" if ($Verbose == 1);
	Leto::GrpAtoms->print_vdw($pdb_name,$path,$Verbose,$Int,$top,$par,$fp,$sp,$ext,"");
	printf (WARN "\n>> Input vdw file inp/$pdb_name\_complex.vdw.inp created from template.\n");
	
	# PROCESSING THE VDW AND ELECTROSTATIC CALCULATION WITH CHARMM
	$begin = time();
	if ($Verbose == 1) {
		warn "\nNow calculating the Van der Waals and electrostatic energies with CHARMM ...\n";
		warn "\tcharmm < inp/$pdb_name\_complex.vdw.inp > out/$pdb_name\_complex.vdw.out\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
	}
	system ("$charmm < inp/$pdb_name\_complex.vdw.inp > out/$pdb_name\_complex.vdw.out");
	$end = time();
	$sec = ($end-$begin) %60;
	$dif = (($end-$begin)-$sec) /60;
	$min = $dif %60;
	($result,$ref_lines) = Leto::GrpAtoms->out_verif("$pdb_name\_complex.vdw",$path,$Verbose);
	foreach $l (0..$#{$ref_lines}) {
		printf (WARN "%s", $ref_lines->[$l]);
	}
	if ($result == 1) {
		warn "$ref_lines->[0]" if ($Verbose == 1);
	}
	else {
		die "\t/!\\$ref_lines->[0]";
	}
	$ref_lines = undef;
	printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min,$sec) if ($Verbose == 1);
	warn "\tSee the log : out/$pdb_name\_complex.vdw.out.\n\n" if ($Verbose == 1);
	printf (WARN "\n>> Calculation of the Van der Waals and electrostatic energies with CHARMM completed in %d seconds.\n\tSee the log : out/%s.vdw.out.\n\n", $sec, $ref_pdb->[-1]);
}


# CREATION OF THE DECOMP_RUN.COM FILE TO PROCEED TO UHBD
if ($CalcUhbd == 1) {
	if (-d "scr/mini") {
		$ext="-mini";
	}
	warn "\nGeneration of the file decomp_run.com from template to process with UHBD on the different chains ...\n" if ($Verbose == 1);
	$boucle=0;
	$grid=0.3; # the grid initial value.
	while ($boucle == 0) { # loop allowing to remake script run.com and calculation if the initial grid were too small.
		Leto::GrpAtoms->print_decomp_run($pdb_name,$Uhbd,$Verbose,$grid,$Int);
		printf (WARN "\n>> UHBD running command file $pdb_name\_complex.decomp_run.com created from template.\n");
	
	# PROCESSING UHBD
		$begin = time();
		foreach $p ("fp","sp","cx") {
			mkpath ("uhbd/elec/$pdb_name\_complex/$pdb_name\_$p/tmp",0,0751);
			copy ("scr/$pdb_name$ext\_$p-cha.crd.his.ter","uhbd/elec/$pdb_name\_complex/$pdb_name\_$p/$pdb_name\_$p-cha.crd");
			copy ("scr/$pdb_name$ext\_$p-rad.crd.his.ter","uhbd/elec/$pdb_name\_complex/$pdb_name\_$p/$pdb_name\_$p-rad.crd");
		}
		if ($Verbose == 1) {
			warn "\nNow executing UHBD to calculate the electrostatic energy ...\n";
			warn "\tcsh inp/$pdb_name\_complex.decomp_run.com > out/$pdb_name\_complex.decomp_run.out\n";
			print "\tPlease wait ";
			sleep(1); print ".";
			sleep(0.5); print ".";
			sleep(0.5); print ".\n";
		}
		system ("csh inp/$pdb_name\_complex.decomp_run.com > out/$pdb_name\_complex.decomp_run.out");
		$end = time();
		$sec = ($end-$begin) %60;
		$dif = (($end-$begin)-$sec) /60;
		$min = $dif %60;
		$dif = ($dif-$min) /60;
		$h = $dif %24;
		($result,$ref_lines) = Leto::GrpAtoms->out_verif("$pdb_name\_complex.decomp_run",$path,$Verbose);
		foreach $l (0..$#{$ref_lines}) {
			printf (WARN "%s", $ref_lines->[$l]);
		}
		if ($result == 1) {
			warn "$ref_lines->[0]" if ($Verbose == 1);
			$boucle=1;
		}
		elsif ($result == 2) {
			$grid=substr($ref_lines->[1],43,5);
			$grid="0".$grid;
			system "rm -r uhbd/elec";
			$boucle=0;
		}
		else {
			die "\t/!\\$ref_lines->[0]";
		}
	
	}
	printf ("Execution of UHBD completed in %d hours %d minutes and %d seconds.\n", $h, $min, $sec) if ($Verbose == 1);
	warn "\tSee the log : out/$pdb_name\_complex.decomp_run.out.\n\n" if ($Verbose == 1);
	printf (WARN "\n>> Execution of UHBD completed in %d hours %d minutes and %d seconds.\n\tSee the log : out/%s\_complex.decomp_run.out.\n\n", $h, $min, $sec, $pdb_name);
}


# CREATION OF THE GBMV FILE TO PROCEED TO GBMV CALCULATION
if ($CalcGbmv == 1) {
	# CREATION OF THE INPUT GBMV FILE
	mkpath ("gbmv",0,0751);

	warn "\nGeneration of the file gbmv.inp from template to process with the GBMV with CHARMM ...\n" if ($Verbose == 1);
	Leto::GrpAtoms->print_gbmv($pdb_name,$path,$Verbose,$Int,$top,$par,$fp,$sp);
		
	printf (WARN "\n>> Input gbmv file inp/$pdb_name\_complex.gbmv.inp created from template.\n");
	# PROCESSING THE GBMV CALCULATION WITH CHARMM
	$begin = time();
	if ($Verbose == 1) {
		warn "\nNow calculating the GBMV energy with CHARMM ...\n";
		warn "\tcharmm < inp/$pdb_name\_complex.gbmv.inp > out/$pdb_name\_complex.gbmv.out\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
	}
	system ("$charmm < inp/$pdb_name\_complex.gbmv.inp > out/$pdb_name\_complex.gbmv.out");
	$end = time();
	$sec = ($end-$begin) %60;
	$dif = (($end-$begin)-$sec) /60;
	$min = $dif %60;
	($result,$ref_lines) = Leto::GrpAtoms->out_verif("$pdb_name\_complex.gbmv",$path,$Verbose);
	foreach $l (0..$#{$ref_lines}) {
		printf (WARN "%s", $ref_lines->[$l]);
	}
	if ($result == 1) {
		warn "$ref_lines->[0]" if ($Verbose == 1);
	}
	else {
		die "\t/!\\$ref_lines->[0]";
	}
	$ref_lines = undef;
	printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min,$sec) if ($Verbose == 1);
	warn "\tSee the log : out/$pdb_name\_complex.gbmv.out.\n\n" if ($Verbose == 1);
	printf (WARN "\n>> Calculation of the elec energy with CHARMM completed in %d minutes and %d seconds.\n\tSee the log : out/%s\_complex.gbmv.out.\n\n", $min, $sec, $pdb_name);
	
	# EXTRACTION OF VALUES AND SUM
	$begin = time();
	if ($Verbose == 1) {
		warn "\nNow calculating the sum energy with a script PERL...\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
	}
	Leto::GrpAtoms->sum_gbmv($pdb_name,$path,$Verbose,$Int,$conf);
	$end = time();
	$sec = ($end-$begin) %60;
	$dif = (($end-$begin)-$sec) /60;
	$min = $dif %60;
	printf ("Execution of PERL script completed in %d minutes and %d seconds.\n", $min,$sec) if ($Verbose == 1);
}


# CALCULATION OF THE FREE ENERGY
if (($CalcSumUhbd == 1) || ($CalcSumGbmv)) {
	# Creation of the complex PDB file name if it doesn't already exist : 
	if ($ref_pdb->[-1] =~ /.+_[a-z]?$/) {
		my ($pdb1,$fc) = split('\_',$ref_pdb->[0]);
		my ($pdb2,$sc) = split('\_',$ref_pdb->[1]);
		push (@{$ref_pdb}, join('_',$pdb1,($fc . $sc)));
	}
	warn "\nCalculation of the free energy as the sum of electrostatic, van der Waals and accessible surface ones ...\n" if ($Verbose == 1);
	Leto::GrpAtoms->energy_sum("uhbd/elec/$pdb_name\_complex/$pdb_name\_cx/tmp/dg-site.dat","scr/$pdb_name.sas.surf","scr/$pdb_name.vdw.inte",$ref_atoms,$pdb_name,@ref);
	warn "\tSee the free energy sum file : out/ener.out.\n\tOne PDB file was created with the free energy for each residue in place of the bfactor field : scr/$pdb_name.ener.pdb.\n\n" if ($Verbose == 1);
	printf (WARN "\n>> Free energy calculated as the sum of electrostatic, van der Waals and accessible surface ones.\n\tSee the free energy sum file : out/ener.out.\n\tOne PDB file was created with the free energy for each residue in place of the bfactor field : scr/$fname.ener.pdb.\n\n");
}


# PLOTTING OF THE RESULTS OF THE FIRST STRUCTURE
if ($FirstPlot == 1) {
	rmtree(['plot/first']);
	mkpath (['plot/first'],0,0751);
	Leto::GrpAtoms->firstplot($ref_pdb);
}


# CREATION OF THE INPUT MINIMIZATION FILE
if ($DynaMini == 1) {
	warn "\nGeneration of the file mini.inp from the template to proceed to the minimization with CHARMM ...\n" if ($Verbose == 1);
	$mini_ext = "his";
	$ref_pdb = Leto::GrpAtoms->print_mini($ref_pdb,$mini_ext,$path,$Verbose,$top,$par);
	printf (WARN "\n>> Input minimization file inp/$ref_pdb->[-1].mini_$mini_ext.inp created from template.\n") if ($Verbose == 1);
	
	# EXECUTION OF THE MINIMIZATION WITH CHARMM
	$begin = time();
	if ($Verbose == 1) {
		warn "\nNow minimizing with CHARMM ...\n";
		warn "\tcharmm < inp/$ref_pdb->[-1].mini_$mini_ext.inp > out/$ref_pdb->[-1].mini_$mini_ext.out\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
	}
	system ("$charmm < inp/$ref_pdb->[-1].mini_$mini_ext.inp > out/$ref_pdb->[-1].mini_$mini_ext.out");
	$end = time();
	$sec = ($end-$begin) %60;
	$dif = (($end-$begin)-$sec) /60;
	$min = $dif %60;
	($result,$ref_lines) = Leto::GrpAtoms->out_verif("$ref_pdb->[-1].mini_$mini_ext",$path,$Verbose);
	foreach $l (0..$#{$ref_lines}) {
		printf (WARN "%s", $ref_lines->[$l]);
	}
	if ($result == 1) {
		warn "$ref_lines->[0]" if ($Verbose == 1);
	}
	else {
		die "\t/!\\$ref_lines->[0]";
	}
	$ref_lines = undef;
	printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min, $sec) if ($Verbose == 1);
	warn "\tSee the log : out/$ref_pdb->[-1].mini_$mini_ext.out.\n\n" if ($Verbose == 1);
	printf (WARN "\n>> Minimization with CHARMM completed in %d minutes and %d seconds.\n\tSee the log : out/%s.mini_$mini_ext.out.\n\n", $min, $sec, $ref_pdb->[-1]) if ($Verbose == 1);
}
else {
	$mini_ext = "his";
}


# CREATION OF THE INPUT DYNAMIC FILES
if ($DynaInp == 1) {
	warn "\nGeneration of the files dyna.inp from the template to proceed to the dynamic with CHARMM ...\n" if ($Verbose == 1);
	($ref_pdb,$restart,$step) = Leto::GrpAtoms->print_dyna($ref_pdb,$mini_ext,$path,$Verbose);
	printf (WARN "\n>> Input dynamic files inp/$ref_pdb->[-1].dyna_prod.inp created from template.\n") if ($Verbose == 1);	
}


# EXECUTION OF THE DYNAMIC WITH CHARMM
if ($DynaExec == 1) {
	$begin = time();
	printf (WARN "Now executing the molecular dynamics with CHARMM ...\n");
	if ($Verbose == 1) {
		warn "\nNow executing the molecular dynamics with CHARMM ...\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
	}
	foreach $d qw(solv heat equil) {
		warn "\tcharmm < inp/$ref_pdb->[-1].dyna_$d.inp > out/$ref_pdb->[-1].dyna_$d.out\n" if ($Verbose == 1);
		printf (WARN "\tcharmm < inp/%s.dyna_%s.inp > out/%s.dyna_%s.out\n", $ref_pdb->[-1], $d, $ref_pdb->[-1], $d);
		system ("$charmm < inp/$ref_pdb->[-1].dyna_$d.inp > out/$ref_pdb->[-1].dyna_$d.out");
		($result,$ref_lines) = Leto::GrpAtoms->out_verif("$ref_pdb->[-1].dyna_$d",$path,$Verbose);
		foreach $l (0..$#{$ref_lines}) {
			printf (WARN "%s", $ref_lines->[$l]);
		}
		if ($result == 1) {
			warn "$ref_lines->[0]" if ($Verbose == 1);
		}
		else {
			die "\t/!\\$ref_lines->[0]";
		}
		$ref_lines = undef;
	}
	for ($i=1;$i<=$restart;$i++) {
		warn "\tcharmm < inp/$ref_pdb->[-1].dyna_prof$i.inp > out/$ref_pdb->[-1].dyna_prod$i.out\n" if ($Verbose == 1);
		printf (WARN "\tcharmm < inp/%s.dyna_%s.inp > out/%s.dyna_%s.out\n", $ref_pdb->[-1], $i, $ref_pdb->[-1], $i);
		system ("$charmm < inp/$ref_pdb->[-1].dyna_prod$i.inp > out/$ref_pdb->[-1].dyna_prod$i.out");
		($result,$ref_lines) = Leto::GrpAtoms->out_verif("$ref_pdb->[-1].dyna_prod$i",$path,$Verbose);
		foreach $l (0..$#{$ref_lines}) {
			printf (WARN "%s", $ref_lines->[$l]);
		}
		if ($result == 1) {
			warn "$ref_lines->[0]" if ($Verbose == 1);
		}
		else {
			die "\t/!\\$ref_lines->[0]";
		}
		$ref_lines = undef;
		warn "\tRestarting ...\n" if (($i>$restart) && ($Verbose == 1));
	}
	$end = time();
	$sec = ($end-$begin) %60;
	$dif = (($end-$begin)-$sec) /60;
	$min = $dif %60;
	$dif = ($dif-$min) /60;
	$h = $dif %24;
	printf ("Execution of the molecular dynamics simulations with CHARMM completed in %d hours %d minutes and %d seconds.\n", $h, $min, $sec) if ($Verbose == 1);
	warn "\tSee the logs : out/$ref_pdb->[-1].dyna_*.out.\n\n" if ($Verbose == 1);
	printf (WARN "\n>> Execution of the molecular dynamics simulations with CHARMM completed in %d hours %d minutes and %d seconds.\n\tSee the logs : out/%s.dyna_*.out.\n\n", $h, $min, $sec, $ref_pdb->[-1]);
}


# Allows to analyze a molecular dynamics which is already made without computing the first steps of the program.
if (($Int == 1) && (!defined $data)) {
	print "\nThe program needs minimum 3 files :\n\tpdb_cx.dyna_prod*.cor\n\tpdb_cx.dyna_desolv.psf\n\tpdb_cx.dyna_desolv.crd or pdb_cx.dyna_solv.crd (if .cor is solvated)\n\tpdb_cx.dyna_solv.psf (if .cor is solvated) \n\n";
	print "Enter the pdb name (name before underscore): ";
	$fname=<STDIN>;
	chomp ($fname);
	print "Enter the name of the first element of the system : ";
	$fc=<STDIN>;
	chomp ($fc);
	print "Enter the name of the second element of the system : ";
	$sc=<STDIN>;
	chomp ($sc);
	
	# Tests if the necessary files are in the directory scr.
	if (!-e glob("scr/$fname\_$fc$sc.dyna_prod*.cor")) {
		if ((!-e "scr/$fname\_$fc$sc.dyna_solv.crd") && (!-e "scr/$fname\_$fc$sc.dyna_desolv.crd")){
			if ((!-e "scr/$fname\_$fc$sc.dyna_solv.psf") && (!-e "scr/$fname\_$fc$sc.dyna_desolv.psf")) {
				die "Warning it misses the following files in directory scr/ :\n $fname\_$fc$sc.dyna_prod*.cor\n $fname\_$fc$sc.dyna_(de)solv.crd\n $fname\_$fc$sc.dyna_(de)solv.psf\n";	
			}
			elsif ((!-e "scr/$fname\_$fc$sc.dyna_desolv.psf") && (-e "scr/$fname\_$fc$sc.dyna_solv.psf")) {
				die "Warning it misses the following files in directory scr/ :\n $fname\_$fc$sc.dyna_prod*.cor\n $fname\_$fc$sc.dyna_(de)solv.crd\n $fname\_$fc$sc.dyna_desolv.psf\n";
			}
			else {
				die "Warning it misses the following files in directory scr/ :\n $fname\_$fc$c.dyna_prod*.cor\n $fname\_$fc$sc.dyna_(de)solv.crd\n";
			}
		}
		elsif ((-e "scr/$fname\_$fc$sc.dyna_solv.crd") && (!-e "scr/$fname\_$fc$sc.dyna_solv.psf")) {
			if (!-e "scr/$fname\_$fc$sc.dyna_desolv.psf") {
				die "Warning it misses the following files in directory scr/ :\n $fname\_$fc$sc.dyna_prod*.cor\n$fname\_$fc$sc.dyna_solv.psf\n$fname\_$fc$sc.dyna_desolv.psf\n"; 
			}
			else {
				die "Warning it misses the following files in directory scr/ :\n $fname\_$fc$sc.dyna_prod*.cor\n$fname\_$fc$sc.dyna_solv.psf\n";
			}
		}
		else {
			if ((!-e "scr/$fname\_$fc$sc.dyna_solv.psf") && (!-e "scr/$fname\_$fc$sc.dyna_desolv.psf")) {
				 die "Warning it misses the following files in directory scr/ :\n $fname\_$fc$sc.dyna_prod*.cor\n $fname\_$fc$sc.dyna_(de)solv.psf\n";
			}
			elsif ((!-e "scr/$fname\_$fc$sc.dyna_desolv.psf") && (-e "scr/$fname\_$fc$sc.dyna_solv.psf")) {
				die "Warning it misses the following files in directory scr/ :\n $fname\_$fc$sc.dyna_prod*.cor\n $fname\_$fc$sc.dyna_desolv.psf\n";
			}
			else {
				 die "Warning it misses the following file in directory scr/ :\n $fname\_$fc$sc.dyna_prod*.cor\n";
			}
		}
	}
	else {
		if ((!-e "scr/$fname\_$fc$sc.dyna_solv.crd") && (!-e "scr/$fname\_$fc$sc.dyna_desolv.crd")) {
			if ((!-e "scr/$fname\_$fc$sc.dyna_solv.psf") && (!-e "scr/$fname\_$fc$sc.dyna_desolv.psf")) {
				die "Warning it misses the following files in directory scr/ :\n $fname\_$fc$sc.dyna_(de)solv.crd\n $fname\_$fc$sc.dyna_(de)solv.psf\n";
			}
			elsif ((!-e "scr/$fname\_$fc$sc.dyna_desolv.psf") && (-e "scr/$fname\_$fc$sc.dyna_solv.psf")) {
				die "Warning it misses the following files in directory scr/ :\n $fname\_$fc$sc.dyna_(de)solv.crd\n $fname\_$fc$sc.dyna_desolv.psf\n";
			}
			else {
				die "Warning it misses the following files in directory scr/ :\n $fname\_$fc$sc.dyna_(de)solv.crd\n";
			}
		}
		elsif ((-e "scr/$fname\_$fc$sc.dyna_solv.crd") && (!-e "scr/$fname\_$fc$sc.dyna_solv.psf")) {
			if (!-e "scr/$fname\_$fc$sc.dyna_desolv.psf") {
				die "Warning it misses the following files in directory scr/ :\n$fname\_$fc$sc.dyna_solv.psf\n$fname\_$fc$sc.dyna_desolv.psf\n"; 
			}
			else {
				die "Warning it misses the following files in directory scr/ :\n$fname\_$fc$sc.dyna_solv.psf\n";
			}	
		}
		else {
			if ((!-e "scr/$fname\_$fc$sc.dyna_solv.psf") && (!-e "scr/$fname\_$fc$sc.dyna_desolv.psf")) {
				die "Warning it misses the following file in directory scr/ :\n $fname\_$fc$sc.dyna_(de)solv.psf\n";
			}
			elsif ((!-e "scr/$fname\_$fc$sc.dyna_desolv.psf") && (-e "scr/$fname\_$fc$sc.dyna_solv.psf")) {
				die "Warning it misses the following files in directory scr/ :\n $fname\_$fc$sc.dyna_desolv.psf\n";
			}
		}
	}
	
	# Allows to know if psf or crd named solv or desolv.
	if ((-e "scr/$fname\_$fc$sc.dyna_solv.psf") && (-e "scr/$fname\_$fc$sc.dyna_solv.crd")) {
		$etat="solv";
	}
	elsif ((-e "scr/$fname\_$fc$sc.dyna_desolv.psf") && (-e "scr/$fname\_$fc$sc.dyna_desolv.crd")) {
		$etat="desolv";
	}
	
	$pdb_files[0]=lc("$fname\_$fc");
	$pdb_files[1]=lc("$fname\_$sc");
	push (@pdb_files,join ('_',$fname,(lc$fc . lc$sc)));
	$ref_pdb=\@pdb_files;
	if (($ConfEner==1) || ($Conf==1) || ($ConfTraj==1)) {
		print "Enter skip used during your dynamic simulation : "; chomp($skip = <STDIN>);
		print "Number of restarts of the dynamics : "; chomp($restart_n = <STDIN>);
		print "Number of coordinate sets in cor file : "; chomp($conf_n = <STDIN>);
	}	
	if (($ConfEner==1) || ($Conf==1) || ($ConfTraj==1) || ($ConfUhbd==1) || ($ConfGbmv==1)) {
		print "\nEnter the name of the topology file (it must be in a lib directory) [default top_all22_prot_na.inp]: ";
		chomp ($top = <STDIN>);
		if ($top eq "") {
			$top="top_all22_prot_na.inp";
		}
		print "Enter the name of the parameter file (it must be in a lib directory) [default par_all22_prot_na.inp]: ";
		chomp ($par = <STDIN>);
		if ($par eq "") {
			$par="par_all22_prot_na.inp";
		}	
	}
	if (($ConfEner==1) || ($Conf==1) || ($ConfTraj==1) || ($ConfUhbd==1) || ($ConfGbmv==1)) {
		print "Enter the definition with charmm keyword of the first part :";
		chomp ($fp=<STDIN>);
		print "Enter the definition with charmm keyword of the second part :";
		chomp ($sp=<STDIN>);
	}
	if (($ConfEner==1) || ($Conf==1)) {
		print "The file tmpsep/enerI.tmp will be open, change the line 56 and put the good lib :\n";
		sleep(3);
		system "vi tmpsep/enerI.tmp";
	}
	if (($ConfTraj==1) || ($Conf==1)) {
		print "The file tmpsep/confI.tmp will be open, change lines 52, 60, 72, 79, 86, 94, 98, 106, 112, 120 and put the good lib :\n";
		sleep(3);
		system "vi tmpsep/confI.tmp";
	}
	if (($ConfUhbd==1) || ($Conf==1)) {
		print "The file tmpsep/sasI.tmp will be open, change lines 49, 53 and put the good lib :\n";
		sleep(3);
		system "vi tmpsep/sasI.tmp";
		print "The file tmpsep/vdwI.tmp will be open, change lines 51, 67, 99 and put the good lib :\n";
		sleep(3);
		system "vi tmpsep/vdwI.tmp";
		print "The file tmpsep/decomp_runI.tmp will be open, put the name of site.dat which must be used :\n";
		sleep(3);
		system "vi tmpsep/decomp_runI.tmp";
	}
	if (($ConfTraj == 1) || ($Conf==1)) {
		print "Number of conformations needed : "; chomp($conf = <STDIN>);		
	}
	#if (($ConfUhbd==1) || ($Conf==1) || ($ConfGbmv==1)) {
		#print "Enter the total number of residue : "; chomp ($res_nb=<STDIN>);
		
		# remove because it isn't very clear for the user (it depends of the vdw.inp)
		#print "Enter the number of residue of the first element of the system : "; chomp ($res1=<STDIN>);
		#print "Enter the number of residue of the second element of the system : "; chomp ($res2=<STDIN>);
	#}
	if (($Conf==1) || ($ConfTraj==1) || ($SumConf==1)) {
		print "Enter the number of chain : "; chomp ($nb_chn=<STDIN>);
		print "Enter the name of chains in the order of appearance in the basic psf separated by a comma : ";chomp ($nom_chn=<STDIN>); 
	}	
		
	if (($Conf==1) || ($ConfTraj==1)) {	
		$fin=0;
		while ($fin==0) {
			print "Enter :  1 if your ligand is form by only one compound\n\t 2 if your ligand is a whole chain\n\t 3 if your ligand is only a part of chain\n\t 4 other\n";chomp ($cas=<STDIN>);
			if ($cas==1) {
				print "Enter the number of residue of the compound (in the original crd): ";chomp ($num_ref=<STDIN>);
				print "Enter the name of the compound (in the original crd): ";chomp ($nam_ref=<STDIN>);$nam_ref=uc $nam_ref;
				$fin=1;
			}
			elsif ($cas==2) {
				print "Enter the name of the chain : \n"; chomp ($chn_ref=<STDIN>);$chn_ref=uc $chn_ref;
				$fin=1;
			}
			elsif ($cas==3) {
				print "Enter the name of the chain : ";chomp ($chn_ref=<STDIN>);$chn_ref=uc $chn_ref;
				print "Enter the number of the first residue of the ligand (in the original crd) : ";chomp ($numfirst_ref=<STDIN>);
				print "Enter the number of the last residue of the ligand (in the original crd) : ";chomp ($numlast_ref=<STDIN>);
				$fin=1;
			}
			elsif ($cas==4) {
				print "Sorry in its current state the program cannot check that the ligand is correctly placed\n";
			}
			else {
				print "You must enter 1,2,3 or 4, start again!\n";
				$fin=0;
			}
		}
	}
	if (($ConfUhbd==1) || ($Conf==1)) {
		print "Enter the name of the site22 that you want use (put it in lib directory) : ";
		$site=<STDIN>;
		chomp ($site);
		print "Enter the size of the grid will be used for the calculation [default 0.4] : ";
		$grid=<STDIN>;
		chomp ($grid);	
		if ($grid eq "") {
			$grid=0.4;
		}
		print "Enter the residue id of the CYS forming disulfide bonds separated by a comma : ";
		$nb_cys=<STDIN>;
		chomp ($nb_cys);
		print "Enter the residue id of the NTER separated with a comma : ";
		$nb_nter=<STDIN>;
		chomp ($nb_nter);
		print "Enter the residue id of the CTER separated with a comma : ";
		$nb_cter=<STDIN>;
		chomp ($nb_cter);
		print "Enter the residue id of the NTNR separated with a comma : ";
		$nb_nter=<STDIN>;
		chomp ($nb_ntnr);
		print "Enter the residue id of the CTNR separated with a comma : ";
		$nb_cter=<STDIN>;
		chomp ($nb_ctnr);
		print "Enter the residue id of the acetylated residu separated with a comma : ";
		$nb_ace=<STDIN>;
		chomp ($nb_ace);
		print "Enter the residue id of the CT3 residu separated with a comma : ";
		$nb_ct3=<STDIN>;
		chomp ($nb_ct3);
		print "Enter the residue id of the CT2 residu separated with a comma : ";
		$nb_ct2=<STDIN>;
		chomp ($nb_ct2);
		print "Enter the residue id of the protonnated asp residu separated with a comma :";
		$nb_app=<STDIN>;
		chomp ($nb_app);
		print "Enter the residue id of the protonnated glu residu separated with a comma :";
		$nb_gpp=<STDIN>;
		chomp ($nb_gpp);
		print "Enter the residue id of the monoanionic phosphoserine separated with a comma : ";
		$nb_sp1=<STDIN>;
		chomp ($nb_sp1);
		print "Enter the residue id of the dianionic phosphoserine separated with a comma : ";
		$nb_sp2=<STDIN>;
		chomp ($nb_sp2);
		print "Enter the nucleic acid id of the 5TER separated with a comma : ";
		$nb_citer=<STDIN>;
		chomp ($nb_citer);
		print "Enter the nucleic acid id of the 3TER acid nucleic separated with a comma : ";
		$nb_tter=<STDIN>;
		chomp ($nb_tter);
		print "Enter the nucleic acid id of the deoxyribose nucleic acid separated with a comma : ";
		$nb_deo=<STDIN>;
		chomp ($nb_deo);
		print "Enter the name of the queue when you want that your UHBD calculation will be submitted : ";
                $queue=<STDIN>;
                chomp ($queue);
		print "Enter your e-mail address to know when your UHBD calculation will be finished : ";
                $mail=<STDIN>;
                chomp ($mail);

	}
}

# Allows to read the informations which are necessary to execute the program from a data file
elsif ((defined $data) && ($Int==1)) {
	open (DATA, "< $data") || die "Unable to open the file $data : $!\n";	
	@dat=<DATA>;
	close(DATA);
	# extracts the necessary variables and removes space
	($fname)=($dat[0]=~/.*:(.*)/);	
	$fname=~s/ *//g;
	($fc)=($dat[1]=~/.*:(.*)/);
	$fc=~s/ *//g;
	($sc)=($dat[2]=~/.*:(.*)/);
	$sc=~s/ *//g;
	($top)=($dat[3]=~/.*:(.*)/);
	$top=~s/ *//g;
	if ($top eq "") {
		$top="top_all22_prot.inp";
	}
	($par)=($dat[4]=~/.*:(.*)/);
	$par=~s/ *//g;
	if ($par eq "") {
		$par="par_all22_prot.inp";
	}
	($skip)=($dat[5]=~/.*:(.*)/);
	$skip=~s/ *//g;
	($restart_n)=($dat[6]=~/.*:(.*)/);
	$restart_n=~s/ *//g;
	($conf_n)=($dat[7]=~/.*:(.*)/);
	$conf_n=~s/ *//g;
	($conf)=($dat[8]=~/.*:(.*)/);
	$conf=~s/ *//g;
	($nb_chn)=($dat[9]=~/.*:(.*)/);
	$nb_chn=~s/ *//g;
	($nom_chn)=($dat[10]=~/.*:(.*)/);
	$nom_chn=~s/ *//g;
	($cas)=($dat[11]=~/.*:(.*)/);
	$cas=~s/ *//g;
	($num_ref)=($dat[16]=~/.*:(.*)/); 
	$num_ref=~s/ *//g;
	($nam_ref)=($dat[17]=~/.*:(.*)/);
	$nam_ref=uc $nam_ref;$nam_ref=~s/ *//g;
	($chn_ref)=($dat[18]=~/.*:(.*)/);
	$chn_ref=uc $chn_ref;$chn_ref=~s/ *//g;
	($numfirst_ref)=($dat[19]=~/.*:(.*)/);
	$numfirst_ref=~s/ *//g;
	($numlast_ref)=($dat[20]=~/.*:(.*)/);
	$numlast=~s/ *//g;
	($fp)=($dat[21]=~/.* bracket :(.*)/);
	($sp)=($dat[22]=~/.* bracket :(.*)/);
	($site)=($dat[24]=~/.*:(.*)/);
	$site=~s/ *//g;
	($grid)=($dat[25]=~/.*:(.*)/);
	$grid=~s/ *//g;
	if ($grid eq "") {
		$grid=0.4;
	}
	($nb_cys)=($dat[26]=~/.*:(.*)/);
	$nb_cys=~s/ *//g;
	($nb_nter)=($dat[27]=~/.*:(.*)/);
	$nb_nter=~s/ *//g;
	($nb_cter)=($dat[28]=~/.*:(.*)/);
	$nb_cter=~s/ *//g;
	($nb_ntnr)=($dat[29]=~/.*:(.*)/);
	$nb_ntnr=~s/ *//g;
	($nb_ctnr)=($dat[30]=~/.*:(.*)/);
	$nb_ctnr=~s/ *//g;
	($nb_ace)=($dat[31]=~/.*:(.*)/);
	$nb_ace=~s/ *//g;
	($nb_ct3)=($dat[32]=~/.*:(.*)/);
	$nb_ct3=~s/ *//g;
	($nb_ct2)=($dat[33]=~/.*:(.*)/);
	$nb_ct2=~s/ *//g;
	($nb_app)=($dat[34]=~/.*:(.*)/);
	$nb_app=~s/ *//g;
	($nb_gpp)=($dat[35]=~/.*:(.*)/);
	$nb_gpp=~s/ *//g;
	($nb_sp1)=($dat[36]=~/.*:(.*)/);
	$nb_sp1=~s/ *//g;
	($nb_sp2)=($dat[37]=~/.*:(.*)/);
	$nb_sp2=~s/ *//g;
	($nb_citer)=($dat[38]=~/.*:(.*)/);
	$nb_citer=~s/ *//g;
	($nb_tter)=($dat[39]=~/.*:(.*)/);
	$nb_tter=~s/ *//g;
	($nb_deo)=($dat[40]=~/.*:(.*)/);
	$nb_deo=~s/ *//g;
	($queue)=($dat[41]=~/.*:(.*)/);
	$queue=~s/ *//g;
	($mail)=($dat[42]=~/.*:(.*)/);
        $mail=~s/ *//g;
        ($nb_nneu)=($dat[43]=~/.*:(.*)/);
        $nb_nneu=~s/ *//g;


	# Allows to know if psf or crd named solv or desolv.
	if ((-e "scr/$fname\_$fc$sc.dyna_solv.psf") && (-e "scr/$fname\_$fc$sc.dyna_solv.crd")) {
		$etat="solv";
	}
	elsif ((-e "scr/$fname\_$fc$sc.dyna_desolv.psf") && (-e "scr/$fname\_$fc$sc.dyna_desolv.crd")) {
		$etat="desolv";
	}
	$pdb_files[0]=lc("$fname\_$fc");
	$pdb_files[1]=lc("$fname\_$sc");
	push (@pdb_files,join ('_',$fname,(lc$fc . lc$sc)));
	$ref_pdb=\@pdb_files;
}

else {
	$restart_n=0;
	$conf_n=0;
	$res_nb=0;
	$etat="solv";
}


# ENERGY CALCULATION FOR EACH CONFORMATION FROM THE MOLECULAR DYNAMICS
if ($ConfEner == 1) {
	if ($Int == 1) {
		# creates the directories inp and out if they don't exist
		mkpath (['inp','out'],0,0751);
		# removes the old directories in the scr
		@old_file= glob ("scr/$ref_pdb->[-1].conf_*");
		rmtree([@old_file]);
	}
	warn "\nGeneration of the file ener.inp from the template to proceed to the energy calculation with CHARMM ...\n" if ($Verbose == 1);
	($ref_pdb,$restart_n,$conf_n) = Leto::GrpAtoms->print_ener($ref_pdb,$path,$Verbose,$Int,$restart_n,$conf_n,$etat,$top,$par,$skip,$fp,$sp);
	printf (WARN "\n>> Input energy calculation file inp/$ref_pdb->[-1].ener.inp created from template.\n") if ($Verbose == 1);
	$begin = time();
	if ($Verbose == 1) {
		warn "\nNow calculating the energy with CHARMM ...\n";
		warn "\t$charmm < inp/$ref_pdb->[-1].ener.inp > out/$ref_pdb->[-1].ener.out\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
	}
	
 	system ("$charmm < inp/$ref_pdb->[-1].ener.inp > out/$ref_pdb->[-1].ener.out");

 	$end = time();
 	$sec = ($end-$begin) %60;
 	$dif = (($end-$begin)-$sec) /60;
 	$min = $dif %60;
 	($result,$ref_lines) = Leto::GrpAtoms->out_verif("$ref_pdb->[-1].ener",$path,$Verbose);
 	foreach $l (0..$#{$ref_lines}) {
 		printf (WARN "%s", $ref_lines->[$l]);
 	}
 	if ($result == 1) {
 		warn "$ref_lines->[0]" if ($Verbose == 1);
 	}
 	else {
 		die "\t/!\\$ref_lines->[0]";
 	}
 	$ref_lines = undef;
 	printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min, $sec) if ($Verbose == 1);
 	warn "\tSee the log : out/$ref_pdb->[-1].ener.out.\n\n" if ($Verbose == 1);
 	printf (WARN "\n>> Execution of CHARMM completed in %d minutes and %d seconds.\n\tSee the log : out/%s.ener.out.\n\n", $min, $sec, $ref_pdb->[-1]);
}


# TRAJECTORIES ANALYSIS
if ($ConfTraj == 1) {
	warn "\nGeneration of the output data file out/$ref_pdb->[-1].ener.dat ...\n" if ($Verbose == 1);
	($ref_pdb,$ref_conf) = Leto::GrpAtoms->traj_analysis($ref_pdb,$conf,$Int);
	printf (WARN "\n>> Data file out/$ref_pdb->[-1].ener.dat created.\n");
	printf (WARN "> Lists in a first part each conformation with its energy.\n");
	printf (WARN "> Classifies into groups the conformations according to their energies and gives the average energy of each group.\n");
# EXTRACTION OF THE SELECTED CONFORMATIONS
	warn "\nGeneration of the file conf.inp from the template to proceed to the extraction of the selected conformations with CHARMM ...\n" if ($Verbose == 1);
	$ref_pdb = Leto::GrpAtoms->print_conf($ref_pdb,$ref_conf,$path,$Verbose,$Int,$conf_n,$restart_n,$etat,$top,$par,$skip,$fp,$sp);
	printf (WARN "\n>> Input dynamic files inp/$ref_pdb->[-1].conf.inp created from template.\n") if ($Verbose == 1);
	$begin = time();
	if ($Verbose == 1) {
		warn "\nNow extracting the selected conformations with CHARMM ...\n";
		warn "\tcharmm < inp/$ref_pdb->[-1].conf.inp > out/$ref_pdb->[-1].conf.out\n";
		print "\tPlease wait ";
		sleep(1); print ".";
		sleep(0.5); print ".";
		sleep(0.5); print ".\n";
	}
	
	system ("$charmm < inp/$ref_pdb->[-1].conf.inp > out/$ref_pdb->[-1].conf.out");
	$end = time();
	$sec = ($end-$begin) %60;
	$dif = (($end-$begin)-$sec) /60;
	$min = $dif %60;
	($result,$ref_lines) = Leto::GrpAtoms->out_verif("$ref_pdb->[-1].conf",$path,$Verbose);
	foreach $l (0..$#{$ref_lines}) {
		printf (WARN "%s", $ref_lines->[$l]);
	}
	if ($result == 1) {
		warn "$ref_lines->[0]" if ($Verbose == 1);
	}
	else {
		die "\t/!\\$ref_lines->[0]";
	}
	$ref_lines = undef;
	printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min, $sec) if ($Verbose == 1);
	warn "\tSee the log : out/$ref_pdb->[-1].conf.out.\n\n" if ($Verbose == 1);
	printf (WARN "\n>> Execution of CHARMM completed in %d minutes and %d seconds.\n\tSee the log : out/%s.conf.out.\n\n", $min, $sec, $ref_pdb->[-1]);

	# Allows to check, before making energy calculation, if in cha and rad.crd files the ligand is positioned in the begin or in the end of the file 
	if ($ref_pdb->[-1] =~ /.+_.?$/) { # Creation of the complex PDB file name if it doesn't already exist
		foreach $f (0..$#{$ref_pdb}) {	
			my ($pdb,$chain) = split('\_',$ref_pdb->[$f]);
			our $fname = $pdb;
			our $complex .= $chain unless ($chain eq "w");
		}
		push (@{$ref_pdb}, join('_',$fname,$complex));
	}
	foreach $b (0..$#{$ref_conf}) { # foreach conformation puts in an array the name.conf* where * is the number of conformation
		$tab[$b] = join('.conf_',$ref_pdb->[-1],$ref_conf->[$b]);
	}
	my $conf_files = \@tab;
	foreach $c (0..$#{$conf_files}) {
		open (FILE, "<scr/$conf_files->[$c]/$conf_files->[$c]-cha.crd.his");
		@line=<FILE>;
		close (FILE);
		$i=$#line;
		chomp($line[$i]);
		chomp($line[6]);
		@last_line=split(' ',$line[$#line]); # puts each field of the crd in an array for the last line of the file
		@first_line=split(' ',$line[6]); # puts each field of the crd in an array for the first line of the file
		if ($cas==1) {
			if ((($first_line[1]==$num_ref) && ($first_line[2] eq $nam_ref)) || (($last_line[1]==$num_ref) && ($last_line[2] eq $nam_ref)))  { # checks if the name and the number of the first or last line corresponds to the ligand
				#print "the ligand is well in the end of file\n";
			}
			else {
				$conf=Leto::GrpAtoms->transfo_crd($cas,$numfirst_ref,$numlast_ref,$num_ref,$nam_ref,$chn_ref,$conf_files->[$c],$ref_pdb); # modifies the crd
				Leto::GrpAtoms->create_crd_for_psf($conf_files->[$c],$cas); # creates crd for each parts of the system (2 proteins, 1 GDP+Mg)
				Leto::GrpAtoms->print_join($cas,$conf_files->[$c],$path,$conf,$Verbose,$top,$par); # allows to creates the inp that make the new psf
				
				$begin = time();
				if ($Verbose == 1) {
					warn "\nNow making psf with CHARMM ...\n";
					warn "\tcharmm < inp/$conf_files->[$c]/$conf_files->[$c].join_crd_ced.inp > out/$conf_files->[$c]/$conf_files->[$c].join_crd_ced.out\n";
					print "\tPlease wait ";
					sleep(1); print ".";
					sleep(0.5); print ".";
					sleep(0.5); print ".\n";
				}
				# Makes the calculation with charmm
				system ("$charmm < inp/$conf_files->[$c]/$conf_files->[$c].join_crd_ced.inp > out/$conf_files->[$c]/$conf_files->[$c].join_crd_ced.out");
				$end = time();
				$sec = ($end-$begin) %60;
				$dif = (($end-$begin)-$sec) /60;
				$min = $dif %60;
				($result,$ref_lines) = Leto::GrpAtoms->out_verif("$conf_files->[$c]/$conf_files->[$c].join_crd_ced",$path,$Verbose);
				foreach $l (0..$#{$ref_lines}) {
					printf (WARN "%s", $ref_lines->[$l]);
				}
				if ($result == 1) {
					warn "$ref_lines->[0]" if ($Verbose == 1);
				}
				else {
					die "\t/!\\$ref_lines->[0]";
				}
				$ref_lines = undef;
				printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min, $sec) if ($Verbose == 1);
				warn "\tSee the log : out/$conf_files->[$c].join_crd_ced.out.\n\n" if ($Verbose == 1);
				printf (WARN "\n>> Execution of CHARMM completed in %d minutes and %d seconds.\n\tSee the log : out/join_crd_ced.out.\n\n", $min, $sec);
			}
			
		}
		
		elsif ($cas==2) {
			if (($first_line[7] eq $chn_ref) || ($last_line[7] eq $chn_ref)) { # checks if the name of the chain of the first or last line corresponds to the ligand
				#print "the ligand is well in the end of file\n";
			}
			else { # if the name of the chain of the first or last line doesn't correspond to the ligand
				$conf=Leto::GrpAtoms->transfo_crd($cas,$numfirst_ref,$numlast_ref,$num_ref,$nam_ref,$chn_ref,$conf_files->[$c],$ref_pdb); # modifies the crd	
				Leto::GrpAtoms->create_crd_for_psf($conf_files->[$c],$cas); # creates crd for each parts of the system (chain B, chain C and chain W)
				Leto::GrpAtoms->print_join($cas,$conf_files->[$c],$path,$conf,$Verbose,$top,$par); # allows to creates the inp that make the new psf
				$begin = time();
				if ($Verbose == 1) {
					warn "\nNow making psf with CHARMM ...\n";
					warn "\tcharmm < inp/$conf_files->[$c]/$conf_files->[$c].join_crd_nico.inp > out/$conf_files->[$c]/$conf_files->[$c].join_crd_nico.out\n";
					print "\tPlease wait ";
					sleep(1); print ".";
					sleep(0.5); print ".";
					sleep(0.5); print ".\n";
				}
				# Makes the calculation with charmm
				system ("$charmm < inp/$conf_files->[$c]/$conf_files->[$c].join_crd_nico.inp > out/$conf_files->[$c]/$conf_files->[$c].join_crd_nico.out");
				$end = time();
				$sec = ($end-$begin) %60;
				$dif = (($end-$begin)-$sec) /60;
				$min = $dif %60;
				($result,$ref_lines) = Leto::GrpAtoms->out_verif("$conf_files->[$c]/$conf_files->[$c].join_crd_nico",$path,$Verbose);
				foreach $l (0..$#{$ref_lines}) {
					printf (WARN "%s", $ref_lines->[$l]);
				}
				if ($result == 1) {
					warn "$ref_lines->[0]" if ($Verbose == 1);
				}
				else {
					die "\t/!\\$ref_lines->[0]";
				}
				$ref_lines = undef;
				printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min, $sec) if ($Verbose == 1);
				warn "\tSee the log : out/$conf_files->[$c].join_crd_nico.out.\n\n" if ($Verbose == 1);
				printf (WARN "\n>> Execution of CHARMM completed in %d minutes and %d seconds.\n\tSee the log : out/join_crd_nico.out.\n\n", $min, $sec);
			}
		}
				
		elsif ($cas==3) {
			if ((($first_line[7] eq $chn_ref) && ($first_line[1] == $numfirst_ref)) || (($last_line[7] eq $chn_ref) && ($last_line[1] == $numlast_ref))) {
				#print "the ligand is well in the end of file\n";
		 	}
			else {
				$conf=Leto::GrpAtoms->transfo_crd($cas,$numfirst_ref,$numlast_ref,$num_ref,$nam_ref,$chn_ref,$conf_files->[$c],$ref_pdb); # modifies the crd
				Leto::GrpAtoms->create_crd_for_psf($conf_files->[$c],$cas); # creates crd for each parts of the system (chain B, chain C and chain W)
				Leto::GrpAtoms->print_join($cas,$conf_files->[$c],$path,$conf,$Verbose,$top,$par); # allows to creates the inp that make the new psf
				$begin = time();
				if ($Verbose == 1) {
					warn "\nNow making psf with CHARMM ...\n";
					warn "\tcharmm < inp/$conf_files->[$c]/$conf_files->[$c].join_crd_ced.inp > out/$conf_files->[$c]/$conf_files->[$c].join_crd_ced.out\n";
					print "\tPlease wait ";
					sleep(1); print ".";
					sleep(0.5); print ".";
					sleep(0.5); print ".\n";
				}
				# Makes the calculation with charmm
				system ("$charmm < inp/$conf_files->[$c]/$conf_files->[$c].join_crd_ced.inp > out/$conf_files->[$c]/$conf_files->[$c].join_crd_ced.out");
				$end = time();
				$sec = ($end-$begin) %60;
				$dif = (($end-$begin)-$sec) /60;
				$min = $dif %60;
				($result,$ref_lines) = Leto::GrpAtoms->out_verif("$conf_files->[$c]/$conf_files->[$c].join_crd_ced",$path,$Verbose);
				foreach $l (0..$#{$ref_lines}) {
					printf (WARN "%s", $ref_lines->[$l]);
				}
				if ($result == 1) {
					warn "$ref_lines->[0]" if ($Verbose == 1);
				}
				else {
					die "\t/!\\$ref_lines->[0]";
				}
				$ref_lines = undef;
				printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min, $sec) if ($Verbose == 1);
				warn "\tSee the log : out/$conf_files->[$c].join_crd_nico.out.\n\n" if ($Verbose == 1);
				printf (WARN "\n>> Execution of CHARMM completed in %d minutes and %d seconds.\n\tSee the log : out/join_crd_ced.out.\n\n", $min, $sec);
			}
		}
		
		else {
			die "Sorry in its current state the program cannot check that the ligand is correctly placed.\n";
		}
		
	}
	if ($cas==4) {
		die "the program does not function yet (functions only in the case of cedric ane nicolas because modifies the crd but it is also necessary to modify the psf)\n";
	}
}

# ALLOWS TO PUT THE GOOD VALUE IN $togsep VARIABLE TO KNOW IF THE USER DOES A CALCULATION WITH ONE SITE PER RESIDUE
# OR ONE SITE FOR BB AND ONE SITE FOR SC
if (($Conf==1) || ($ConfSas==1) || ($ConfVdw==1) || ($ConfUhbd==1) || ($SumConf == 1) || ($average == 1) || ($stdev == 1) || ($average_nonpond == 1) || ($ConfPlot == 1) ) {
	if ($Separated == 1) {
		$togsep="sep";
	}
	elsif ($Together == 1) {
		$togsep="tog";
	}
	print $togsep;
}		



# CALCULATIONS OF SAS ON THE SELECTED CONFORMATIONS
if ($ConfSas == 1) {
	if ($ref_pdb->[-1] =~ /.+_.?$/) {
		foreach $f (0..$#{$ref_pdb}) {	
			my ($pdb,$chain) = split('\_',$ref_pdb->[$f]);
			our $fname = $pdb;
			our $complex .= $chain unless ($chain eq "w");
		}
		push (@{$ref_pdb}, join('_',$fname,$complex));
	}
	my @temp = glob("scr/$ref_pdb->[-1].conf_*");
	foreach $t (0..$#temp) { # foreach conformation puts in a array the name.conf* where * is the number of conformation
		my ($rep,$name,$conf) = split(/\/|\./,$temp[$t]);
		$tab[$t] = join('.',$name,$conf);
	}

	my $conf_files = \@tab;
	foreach $c (0..$#{$conf_files}) {
		my $ref_pdb_conf = $ref_pdb;
		push (@{$ref_pdb_conf},"$conf_files->[$c]");
	
		# CREATION OF THE INPUT SAS FILE
		warn "\nGeneration of the file sas.inp from template to process with the sas calculation with CHARMM ...\n" if ($Verbose == 1);
		Leto::GrpAtoms->print_sas($ref_pdb_conf,$path,$Verbose,$Int,$top,$par,$fp,$sp,"",$togsep);
		printf (WARN "\n>> Input sas file inp/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].sas.inp created from template.\n");
		# PROCESSING THE SAS CALCULATION WITH CHARMM
		$begin = time();
		if ($Verbose == 1) {
			warn "\nNow calculating the solvent accessible surface energy with CHARMM...\n";
			warn "\tcharmm < inp/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].sas.inp > out/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].sas.out\n";
			print "\tPlease wait ";
			sleep(1); print ".";
			sleep(0.5); print ".";
			sleep(0.5); print ".\n";
		}
		system ("$charmm < inp/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].sas.inp > out/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].sas.out");
		$end = time();
		$sec = ($end-$begin) %60;
		$dif = (($end-$begin)-$sec) /60;
		$min = $dif %60;
		($result,$ref_lines) = Leto::GrpAtoms->out_verif("$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].sas",$path,$Verbose);
		foreach $l (0..$#{$ref_lines}) {
			printf (WARN "%s", $ref_lines->[$l]);
		}
		if ($result == 1) {
			warn "$ref_lines->[0]" if ($Verbose == 1);
		}
		else {
			die "\t/!\\$ref_lines->[0]";
		}
		$ref_lines = undef;
		printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min,$sec) if ($Verbose == 1);
		warn "\tSee the log : out/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].sas.out.\n\n" if ($Verbose == 1);
		printf (WARN "\n>> Calculation of the solvent accessible surface energy with CHARMM completed in %d minutes and %d seconds.\n\tSee the log : out/%s.sas.out.\n\n", $min, $sec, $ref_pdb_conf->[-1]);
	}
}	


# CALCULATIONS OF VDW ON THE SELECTED CONFORMATIONS
if ($ConfVdw == 1) {
	if ($ref_pdb->[-1] =~ /.+_.?$/) {
		foreach $f (0..$#{$ref_pdb}) {	
			my ($pdb,$chain) = split('\_',$ref_pdb->[$f]);
			our $fname = $pdb;
			our $complex .= $chain unless ($chain eq "w");
		}
		push (@{$ref_pdb}, join('_',$fname,$complex));
	}
	my @temp = glob("scr/$ref_pdb->[-1].conf_*");
	foreach $t (0..$#temp) { # foreach conformation puts in a array the name.conf* where * is the number of conformation
		my ($rep,$name,$conf) = split(/\/|\./,$temp[$t]);
		$tab[$t] = join('.',$name,$conf);
	}

	my $conf_files = \@tab;
	foreach $c (0..$#{$conf_files}) {
		my $ref_pdb_conf = $ref_pdb;
		push (@{$ref_pdb_conf},"$conf_files->[$c]");
	
		# CREATION OF THE INPUT VDW FILE
		warn "\nGeneration of the file vdw.inp from template to process with the Van der Waals and electrostatic energies calculation with CHARMM ...\n" if ($Verbose == 1);
		Leto::GrpAtoms->print_vdw($ref_pdb_conf,$path,$Verbose,$Int,$top,$par,$fp,$sp,"",$togsep);
		printf (WARN "\n>> Input vdw file inp/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].vdw.inp created from template.\n");
	
		# PROCESSING THE VDW AND ELECTROSTATIC CALCULATION WITH CHARMM
		$begin = time();
		if ($Verbose == 1) {
			warn "\nNow calculating the Van der Waals and electrostatic energies with CHARMM ...\n";
			warn "\tcharmm < inp/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].vdw.inp > out/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].vdw.out\n";
			print "\tPlease wait ";
			sleep(1); print ".";
			sleep(0.5); print ".";
			sleep(0.5); print ".\n";
		}
		system ("$charmm < inp/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].vdw.inp > out/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].vdw.out");
		$end = time();
		$sec = ($end-$begin) %60;
		$dif = (($end-$begin)-$sec) /60;
		$min = $dif %60;
		($result,$ref_lines) = Leto::GrpAtoms->out_verif("$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].vdw",$path,$Verbose);
		foreach $l (0..$#{$ref_lines}) {
			printf (WARN "%s", $ref_lines->[$l]);
		}
		if ($result == 1) {
			warn "$ref_lines->[0]" if ($Verbose == 1);
		}
		else {
			die "\t/!\\$ref_lines->[0]";
		}
		$ref_lines = undef;
		printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n",$min,$sec) if ($Verbose == 1);
		warn "\tSee the log : out/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].vdw.out.\n\n" if ($Verbose == 1);
		printf (WARN "\n>> Calculation of the Van der Waals and electrostatic energies with CHARMM completed in %d seconds.\n\tSee the log : out/%s.vdw.out.\n\n", $sec, $ref_pdb_conf->[-1]);
		
	}
}


# CALCULATIONS OF PB ELEC ON THE SELECTED CONFORMATIONS USING ONE TITRABLE SITE PER RESIDUE
if ($ConfUhbd == 1) {
	if ($ref_pdb->[-1] =~ /.+_.?$/) {
		foreach $f (0..$#{$ref_pdb}) {	
			my ($pdb,$chain) = split('\_',$ref_pdb->[$f]);
			our $fname = $pdb;
			our $complex .= $chain unless ($chain eq "w");
		}
		push (@{$ref_pdb}, join('_',$fname,$complex));
	}
	my @temp = glob("scr/$ref_pdb->[-1].conf_*");
	foreach $t (0..$#temp) { # foreach conformation puts in a array the name.conf* where * is the number of conformation
		my ($rep,$name,$conf) = split(/\/|\./,$temp[$t]);
		$tab[$t] = join('.',$name,$conf);
	}

	my $conf_files = \@tab;
	foreach $c (0..$#{$conf_files}) {
		my $ref_pdb_conf = $ref_pdb;
		push (@{$ref_pdb_conf},"$conf_files->[$c]");
		
		# CREATION OF THE DECOMP_RUN.COM FILE TO PROCEED TO UHBD
		warn "\nGeneration of the file decomp_run.com from template to process with UHBD on the different chains ...\n" if ($Verbose == 1);
		#$grid=0.40;
		Leto::GrpAtoms->print_decomp_run($ref_pdb_conf,$Uhbd,$Verbose,$grid,$Int,$site,$path);
		printf (WARN "\n>> UHBD running command file $ref_pdb_conf->[-1].decomp_run.com created from template.\n");
		my ($name,$conf) = ($ref_pdb_conf->[-1]=~m/(.*)\.(.*)/);
		pop (@{$ref_pdb});
		foreach $p (0..$#{$ref_pdb}) {
			if ($ref_pdb->[$p] !~ /._w$/) {
				mkpath ("uhbd/elec/$conf_files->[$c]/$ref_pdb->[$p].$conf/tmp",0,0751);
				copy ("scr/$conf_files->[$c]/$ref_pdb->[$p].$conf-cha.crd.his","uhbd/elec/$conf_files->[$c]/$ref_pdb->[$p].$conf/$ref_pdb->[$p].$conf-cha.crd");
				copy ("scr/$conf_files->[$c]/$ref_pdb->[$p].$conf-rad.crd.his","uhbd/elec/$conf_files->[$c]/$ref_pdb->[$p].$conf/$ref_pdb->[$p].$conf-rad.crd");
				warn "\nCorrection of the CRD files to comply with the UHBD standard ...\n" if ($Verbose == 1);
				if ($Int == 0) {
					Leto::GrpAtoms->crd2uhbd($ref_atoms,$ref_last,"uhbd/elec/$conf_files->[$c]/$ref_pdb->[$p].$conf/*-{cha,rad}.crd");
				}
				else {
					Leto::GrpAtoms->crd2uhbd_interactif("uhbd/elec/$conf_files->[$c]/$ref_pdb->[$p].$conf/*-{cha,rad}.crd",$nb_cys,$nb_nter,$nb_cter,$nb_ntnr,$nb_ctnr,$nb_ace,$nb_ct3,$nb_ct2,$nb_app,$nb_gpp,$nb_sp1,$nb_sp2,$nb_citer,$nb_tter,$nb_deo,$nb_nneu);
				}
				printf (WARN "\n>> CRD files corrected to be used in UHBD.\n");
				printf (WARN "> Backbone atoms of the first and the last residues of each chain renamed in NTER or CTER.\n");
				printf (WARN "> Atoms of the protonnated ASP and GLU residues renamed in APP and GPP, respectively.\n");
				printf (WARN "> Atoms of the cysteins involved in disulfide bonds renamed in CSS.\n");
			}
		}
	}
}


# CALCULATIONS OF GB ELEC ON THE SELECTED CONFORMATIONS
if ($ConfGbmv == 1) {
	if ($ref_pdb->[-1] =~ /.+_.?$/) {
		foreach $f (0..$#{$ref_pdb}) {	
			my ($pdb,$chain) = split('\_',$ref_pdb->[$f]);
			our $fname = $pdb;
			our $complex .= $chain unless ($chain eq "w");
		}
		push (@{$ref_pdb}, join('_',$fname,$complex));
	}
	my @temp = glob("scr/$ref_pdb->[-1].conf_*");
	foreach $t (0..$#temp) { # foreach conformation puts in a array the name.conf* where * is the number of conformation
		my ($rep,$name,$conf) = split(/\/|\./,$temp[$t]);
		$tab[$t] = join('.',$name,$conf);
	}
	
	my $conf_files = \@tab;
	foreach $c (0..$#{$conf_files}) {
		my $ref_pdb_conf = $ref_pdb;
		push (@{$ref_pdb_conf},"$conf_files->[$c]");
		
		# CREATION OF THE INPUT GBMV FILE
		foreach $p (0..$#{$ref_pdb}) {
			if ($ref_pdb->[$p] !~ /._w$/) {
				mkpath ("gbmv/$conf_files->[$c]",0,0751);
			}
		}
		warn "\nGeneration of the file gbmv.inp from template to process with the GBMV with CHARMM ...\n" if ($Verbose == 1);
		
		Leto::GrpAtoms->print_gbmv($ref_pdb_conf,$path,$Verbose,$Int,$top,$par,$fp,$sp);
		
		printf (WARN "\n>> Input gbmv file inp/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].gbmv.inp created from template.\n");
		my ($name,$conf) = ($ref_pdb_conf->[-1]=~m/(.*)\.(.*)/);
		
		
		# PROCESSING THE GBMV CALCULATION WITH CHARMM
		$begin = time();
		if ($Verbose == 1) {
			warn "\nNow calculating the GBMV energy with CHARMM ...\n";
			warn "\tcharmm < inp/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].gbmv.inp > out/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].gbmv.out\n";
			print "\tPlease wait ";
			sleep(1); print ".";
			sleep(0.5); print ".";
			sleep(0.5); print ".\n";
		}
		system ("$charmm < inp/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].gbmv.inp > out/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].gbmv.out");
		$end = time();
		$sec = ($end-$begin) %60;
		$dif = (($end-$begin)-$sec) /60;
		$min = $dif %60;
		($result,$ref_lines) = Leto::GrpAtoms->out_verif("$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].gbmv",$path,$Verbose);
		foreach $l (0..$#{$ref_lines}) {
			printf (WARN "%s", $ref_lines->[$l]);
		}
		if ($result == 1) {
			warn "$ref_lines->[0]" if ($Verbose == 1);
		}
		else {
			die "\t/!\\$ref_lines->[0]";
		}
		$ref_lines = undef;
		printf ("Execution of CHARMM completed in %d minutes and %d seconds.\n", $min,$sec) if ($Verbose == 1);
		warn "\tSee the log : out/$ref_pdb_conf->[-1]/$ref_pdb_conf->[-1].gbmv.out.\n\n" if ($Verbose == 1);
		printf (WARN "\n>> Calculation of the elec energy with CHARMM completed in %d minutes and %d seconds.\n\tSee the log : out/%s.gbmv.out.\n\n", $min, $sec, $ref_pdb_conf->[-1]);
	
		# EXTRACTION OF VALUES AND SUM
		$begin = time();
		if ($Verbose == 1) {
			warn "\nNow calculating the sum energy with a script PERL...\n";
			print "\tPlease wait ";
			sleep(1); print ".";
			sleep(0.5); print ".";
			sleep(0.5); print ".\n";
		}
		$ref_pdb_conf = Leto::GrpAtoms->sum_gbmv($ref_pdb_conf,$path,$Verbose,$Int,$conf);
		$end = time();
		$sec = ($end-$begin) %60;
		$dif = (($end-$begin)-$sec) /60;
		$min = $dif %60;
		printf ("Execution of PERL script completed in %d minutes and %d seconds.\n", $min,$sec) if ($Verbose == 1);
	}
}


# Creates the files to submit uhbd jobs on hpc
if ($Run ==1) {
	# Opens the file that contains the number of conformation
	open (CONF, "< out/$fname\_$fc$sc.distrib.ener.dat");
	chomp(@line_conf=<CONF>);
	close(CONF);
	# The number of conformation is contained in the first column of the file out/name_cx.distrib.ener.dat
	foreach $i(1..$#line_conf) {
		@line=split('\s',$line_conf[$i]);
		push (@conf,$line[0]);
	}
	# Creates one file per conformation
        foreach $a(0..$#conf) {
                open (RUN, "> $conf[$a].run");
               printf (RUN "#!/bin/bash -f\n");
               printf (RUN "#PBS -A cbu\n");
	       printf (RUN "#PBS -q idle\n");
               printf (RUN "#PBS -N \"decomp_$conf[$a]\" \n");
               printf (RUN "#PBS -l walltime=12:00:00\n");
               printf (RUN "#PBS -l nodes=1:ppn=2 \n");
               printf (RUN "#PBS -l pmem=800mb \n");
               printf (RUN "date\n\n");
               printf (RUN "cd $path2\n\n");
               printf (RUN "./inp/$fname\_$fc$sc.conf_$conf[$a]/$fname\_$fc$sc.conf_$conf[$a].decomp_run.com  > out/$fname\_$fc$sc.conf_$conf[$a]/$fname\_$fc$sc.conf_$conf[$a].decomp_run.out\n");
        }
}


# CALCULATION OF THE FREE ENERGY
if ($SumConf == 1) {
	my $atom,$atomnum,$atomname,$resname,$resnum,$x,$y,$z,$occ,$bfactor,$chaine,$li;
	my @lines;
	if ($ref_pdb->[-1] =~ /.+_.?$/) {
		foreach $f (0..$#{$ref_pdb}) {  
 			my ($pdb,$chain) = split('\_',$ref_pdb->[$f]);
 			our $fname = $pdb;
 			our $complex .= $chain unless ($chain eq "w");
 		}
		push (@{$ref_pdb}, $fname."_".$complex);		
	}

	my @temp = glob("scr/$ref_pdb->[-1].conf_*/$ref_pdb->[-1].conf_*.sas.surf");
	foreach $t (0..$#temp) {
		my ($rep,$name,$conf,$ext) = split(/\/|\./,$temp[$t]);
		$tab[$t] = join('.',$name,$conf);
	}
	my $conf_files = \@tab;
	
	foreach $c (0..$#{$conf_files}) {
	
	#Coordinates of each conformation are copied in a new objet new_atoms
	#each new_atoms put in an array which contains adress of each new_atoms.
	#fonction read-pdb can't be used to read pdb file created by the programm
	#and coordinates can't be modified easily in the ref_atoms created when the original
	#pdb file is read because they don't contain the same number of atoms (hydrogens)
		open (FILIN, "< scr/$conf_files->[$c]/$conf_files->[$c].pdb") || die "Unable to open the file $conf_files->[$c]/$conf_files->[$c].pdb : $!\n";
		my @lines = <FILIN>;
		close (FILIN);
		$li=0; 
		undef @new_ref_atoms;
	#coordinates of the pdb file of each conformation are cleared before a new conformation  
		foreach $l (0..$#lines) {
			if ($lines[$l] =~ /^ATOM/) {
				($atom,$atomnum,$atomname,$resname,$resnum,$x,$y,$z,$occ,$bfactor,$chaine) = split(' ',$lines[$l]);
                  		my $new_atom = Leto::Atom->new();
				$new_atom->{RECORD_NAME} = $atom;
				$new_atom->{ATOM_NUMBER} = $atomnum;
				$new_atom->{ATOM_NAME} = $atomname;
				$new_atom->{RESIDUE_NAME} = $resname;
				$new_atom->{RESIDUE_NUMBER} = $resnum;  		
				$new_atom->{X} = $x;
				$new_atom->{Y} = $y;
				$new_atom->{Z} = $z;
				$new_atom->{CHAIN} = $chaine;			
				$li++;
				push(@new_ref_atoms,$new_atom); 		
			}		  		
               }
		  
		# READING OF THE FILE
		warn "\nCalculation of the free energy as the sum of electrostatic, van der Waals and accessible surface ones ...\n" if ($Verbose == 1);
		if (-e "uhbd/elec/$conf_files->[$c]/$conf_files->[$c]/tmp/dg-site.dat") {
			Leto::GrpAtoms->energy_sum_conf("uhbd","uhbd/elec/$conf_files->[$c]/$conf_files->[$c]/tmp/dg-site.dat","scr/$conf_files->[$c]/$conf_files->[$c].sas.surf","scr/$conf_files->[$c]/$conf_files->[$c].vdw.inte",\@new_ref_atoms,$ref_pdb,$conf_files->[$c],$nb_chn,$nom_chn,$togsep,@ref);
		}
		elsif (-e "gbmv/$conf_files->[$c]/elec.crd") {
			Leto::GrpAtoms->energy_sum_conf("gbmv","gbmv/$conf_files->[$c]/elec.crd","scr/$conf_files->[$c]/$conf_files->[$c].sas.surf","scr/$conf_files->[$c]/$conf_files->[$c].vdw.inte",\@new_ref_atoms,$ref_pdb,$conf_files->[$c],$nb_chn,$nom_chn,$togsep,@ref);
		}
		else {
			die "file uhbd/elec/$conf_files->[$c]/$conf_files->[$c]/tmp/dg-site.dat or gbmv/$conf_files->[$c]/elec.crd don't exist\n";
		}
		warn "\tSee the free energy sum file : out/ener.out.\n\tOne PDB file was created with the free energy for each residue in place of the bfactor field : scr/$fname.ener.pdb.\n\n" if ($Verbose == 1);
		printf (WARN "\n>> Free energy calculated as the sum of electrostatic, van der Waals and accessible surface ones.\n\tSee the free energy sum file : out/ener.out.\n\tOne PDB file was created with the free energy for each residue in place of the bfactor field : scr/$fname.ener.pdb.\n\n");
	}
}


# CALCULATION OF THE WEIGHTED (OR NO WEIGHTED) AVERAGE AND STANDARD DEVIATION
if (($average == 1) || ($stdev == 1) || ($average_nonpond == 1)) {
	$ave=0;
	$std=0;
	$ave_nonpond=0;
	if ($average == 1) {
		$ave=1;
	}
	if ($stdev == 1) {
		$std=1;
	}
	if ($average_nonpond == 1) {
		$ave_nonpond=1;
	}
	Leto::GrpAtoms->average_stdev($ref_pdb,$ave,$std,$ave_nonpond,$togsep);
}


# PLOTTING OF THE RESULTS OF THE AVERAGE OF SEVERAL CONFORMATIONS
if ($ConfPlot == 1) {
	rmtree(['plot/average']);
	mkpath (['plot/average'],0,0751);
	rmtree(['plot/efficiency-average']);
	mkpath (['plot/efficiency-average'],0,0751);
	Leto::GrpAtoms->confplot($ref_pdb,$togsep);
}


# ADDITIONAL FUNCTION
if ($supplement ==1) {
	Leto::GrpAtoms->supplement($ref_pdb);
}


# END OF THE PROGRAMM 
$stop = time();
$sec = ($stop-$start) %60;
$dif = (($stop-$start)-$sec) /60;
$min = $dif %60;
$dif = ($dif-$min) /60;
$h = $dif %24;
printf ("\nProcessing of the file %s successfully completed in %d hours %d minutes and %d seconds. (see out/%s.out for more informations)\n\n", $filin, $h, $min, $sec, $fname) if ($Verbose == 1);
printf (WARN "\n>> Processing of the file %s successfully completed in %d hours %d minutes and %d seconds.\n", $filin, $h, $min, $sec);
printf (WARN "\nFor more information about the way this script works, please refer to its documentation and its appending module Leto.pm.\n\n");
close (WARN);
