
# Reading ISIS base (Reaccs) input data
# data format : see JCICS 32(3) 1992 244

require HashObject ;

# HISTORY :
#	BB 28-3-95
# 	Reading the "$MOL" kw not in MDL_molecule class but in
#    MDL_reaction_component.

# Change BB 8-4-95
# reading/storing the header block for molecules
#  e.g. -ISIS-  01179509292D

# Change BB Dec 2 1998
# correct header block
# polar hydrogen removal

$[ = 1 ;

#**************************************************************************#
package MDL_atom ;
#**************************************************************************#
	@ISA = qw( HashObject );

	@fields_names = ('x', 'y', 'z', 'symbol', 'mass_difference', 'charge', 'atom_stereo_parity', 'hydrogens_counts_plus_1', 'stereo_care_box', 'valence', 'HO_designator', 'reaction_component_type', 'reaction_component_number', 'atom_atom_mapping_number', 'inversion_retention_flag', 'exact_change_flag') ;

	eval( ${HashObject'new_SubClass_String} ) ;

#--------------------------------------------------------------------------#

	sub readFromCurrentLine {
			my $class = shift ;
			my $self = new($class) ;
			$self || die 'Assertion failed : Undefined value' ;
			my @fields = /(..........)(..........)(..........).(...)(..)(...)(...)?(...)?(...)?(...)?(...)?(...)?(...)?(...)?(...)?(...)?/ ;
	
		#print STDERR "Fields .. = $#fields\n";
		if( $#fields < 6 ) {
				
			die( "Atom line parsing error, not enough fields: $#fields fields\n $_\n " );
		}
		
		
		foreach $name (@fields_names)	{
			$self->{$name} = shift @fields ;
		}
		
		#remove blanks from symbol name
		$self->{'symbol'} =~ s/\s+//g ;
		$self; 		

	}

#--------------------------------------------------------------------------#
		
	sub	writeLine	{	
                my $self = shift;
				my @fields = () ;
				
				foreach $name (@fields_names)	{
					push(@fields, $self->{$name}) ;
				}

				printf "%10.4f%10.4f%10.4f %-3s%2d%3d%3d%3d%3d%3d%3d%3d%3d%3d%3d%3d\n", @fields ;
	}


#**************************************************************************#
package MDL_bond ;
#**************************************************************************#

	@ISA = qw( HashObject );

	@fields_names = ('first_atom', 'second_atom', 'type', 'stereo', 'not_used', 'topology', 'reacting_center') ;

	eval( ${HashObject'new_SubClass_String} ) ;
#--------------------------------------------------------------------------#
	
	sub readFromCurrentLine {
			my $class = shift ;
			my $self = new($class) ;

# modified for Quanta : the MDL exported has one number less		
#		print STDERR $#fields, "\n" ;
		my @fields = /^([0-9 ][0-9 ][0-9])([0-9 ][0-9 ][0-9])([0-9 ][0-9 ][0-9])([0-9 ][0-9 ][0-9])([0-9 ][0-9 ][0-9])?([0-9 ][0-9 ][0-9])?([0-9 ][0-9 ][0-9])?/ ;

		if($#fields !=  $#fields_names ) {
			die( "Bond line parsing error : \n $_\n " );
		}
		
		foreach $name (@fields_names)	{
			my $value = shift @fields ;
			if(!defined $value) {$value = 0} ;
			$self->{$name} = $value ;
		}
		$self; 		
	}

#--------------------------------------------------------------------------#
		
sub	writeLine	{	
	my $self = shift;
	my @fields = () ;
	
	my $saved;
	
	
#	# smallest index must coming first ?
#	if( $self->{'first_atom'} > $self->{'second_atom'} )	{
#		$saved = $self->{'first_atom'} ;
#		$self->{'first_atom'} = $self->{'second_atom'} ;
#		$self->{'second_atom'} = $saved ;				
#	}
	
	# Some checking
	if( $self->{'first_atom'} == $self->{'second_atom'} )	{
		die("There is bond between the atom $self->{'first_atom'} and itself") ;			
	}
	
	
	foreach $name (@fields_names)	{
		push(@fields, $self->{$name}) ;
	}

	printf( "%3d%3d%3d%3d%3d%3d%3d\n", @fields) || die "Write error" ;
}
#--------------------------------------------------------------------------#

sub connectsAtomIndexes {
	my $self = shift;
	my $a1 = shift ;
	my $a2 = shift ;

	$a1 != $a2 || die "Assertion failed" ;
	
	if( $a1 == $self->{'first_atom'} && $a2 == $self->{'second_atom'} ) {
		return 1 ;
	}

	if( $a2 == $self->{'first_atom'} && $a1 == $self->{'second_atom'} ) {
		return 1 ;
	}

	return 0 ;

}
#**************************************************************************#
package MDL_atom_list_line ;
#**************************************************************************#
	@ISA = qw( HashObject );

	@fields_names = ('atom_number', 'k', 'n_entries', 'entries', ) ;

	eval( ${HashObject'new_SubClass_String} ) ;

#--------------------------------------------------------------------------#

sub readFromCurrentLine 
{
	my $class = shift ;
	my $self = new($class) ;
	$self || die 'Assertion failed : Undefined value' ;
# 27 F    5  17   9  35  53  85
#   4 F    2   7   8
#Changed Sept 15 th 1999 : more flexible, cfr test_3d_query
	my @fields = /^\s*(\d{1,3})\s+([TF])\s{4}([0-5])/ ;
	@fields = (@fields , split " ", $'); #Split like awk: p.221
	#print STDERR "Fields  = $#fields\n", join(":", @fields), "\n";
	if( $#fields < 2 ) {
		 die( "Atom list line parsing error, not enough fields: $#fields fields\n $_\n " );
	}

	

	foreach $name ('atom_number', 'k', 'n_entries')	{
		$self->{$name} = shift @fields ;
	}
	$self->{'entries'} = \@fields ;

	($#fields == $self->{'n_entries'}) || die "Not enough entries:\n$_\n";
	 
	$self; 		

 }

#--------------------------------------------------------------------------#
sub	writeLine
{	
	my $self = shift;
	my @fields = () ;

	foreach $name (@fields_names)	{
		push(@fields, $self->{$name}) ;
	}
	
	printf( "%3d %s    %1d", @fields) || die "Write error" ;
	
	foreach $number ( @{$self->{'entries'}})
	{
		printf( "%4d", $number) || die "Write error" ;
	} 
	
	print( "\n") || die "Write error" ;

				
}


#**************************************************************************#
package MDL_non_parsed_molecule ; 
#useful for quick sdf parsing 
#The connection table is not parsed but stored in a string
#**************************************************************************#

	@ISA = qw( HashObject );

	@fields_names = (

	#new :3 lines of the header block
	'line1', #molecule name
	'line2', # blank or a complex line
	'line3', # comment
	
	'number_of_atoms',
	'number_of_bonds',
	'atoms',
	'bonds',
	'atom_list_block',
	'stext_block',
	
	'count_line',
	#'header_block', #obsolete
	 'properties_block',
	 'data',
	 'id',
	 'original_count_line',

	#For the Count line:
	 'number_of_atoms_list',
	 'obsolete',
	 'chiral_flag',
	 'number_of_stext_entries',
	 'number_of_reaction_components',
	 'number_of_reactants',
	 'number_of_products',
	 'number_of_intermediates',
	 'number_of_add_prop_lines',
	 'ctab_version',
	 
	
	#To store the molecule in one string for quick printing/parsing
	'mol_string'
	 )
	  ;

	# obsolete
	# $header line is a line like   "-ISIS-  01069517512D " which appears
	# at the beginning of the MDL mol file

	# $header block #obsolete
	
	# @properties_block is an array. Each element contains
	# a line of the MDL mol file which comes after the bond block
	# to the "M  END" line

	# %data is an associative array containing MW, CAS, SYMBOL and MOLREGNO

	# $id is normally = $self->{'data'}{'CAS'} when this latter is defined
	# the link is made in the function 
	#              MDL_Reaction'link_molecules_with_data_blocks

	eval( ${HashObject'new_SubClass_String} ) ;
#--------------------------------------------------------------------------#
#From page 2-8 from the new PDF manual
sub version_stamp_is_present {
	my $self = shift ;
		
	return $self->{'ctab_version'} =~ /^\s+v/i;
}
#--------------------------------------------------------------------------#
sub readFromInput {
		my $class = shift ;
		my $self = new($class) ;
		my @atoms = () ;
		my @bonds = () ;
		my $d = "[0-9]+" ;
		my $j  ;
		
# Change BB 28-3-95
# put in MDL_reaction_components
		
#		while( <STDIN> )	{
#			if( /^\$MOL/ )	{ last }
#		}

	#read the header block
		foreach $line (qw(line1 line2 line3))
		{
			$_ = <STDIN> ;
			defined $_ || return undef ;
			chomp ;
			$self->{$line} = $_ ;
		}
	
	#read the count line
		#  6  5  0  0  0  0  0  0  0  0  1 
		#    or 
		# 10  9  0  0  0  0  0  0  0  0  1 V2000

		$_ = <STDIN> ;
		defined $_ || return undef ;		
#
# Changed to be compatible with mdl files
# exported from Quanta
#			if(! /(..\d)(..\d)(\s+$d\s+$d\s+$d\s+$d)/o ) 
		if( ! /(..\d)(..\d)(.*)/o ) #NCI files have only two fields
		{
			print STDERR "Current input line is \n", $_, "\n" ;
			
			die "File format error" ;
			return undef ;
		}
		$self->{'number_of_atoms'} = $1 ;
		$self->{'number_of_bonds'} = $2 ;
		$self->{'count_line'}	= $3 ;
		 
		{
		my @fields = /..\d..\d(...)(...)(...)(...)(...)(...)(...)(...)(...)(......)/;
#		print STDERR "Fields = ",$#fields, "\n" ;
		if( $#fields >= 0)
		{
			$self->{'number_of_atoms_list'} = shift @fields ;
			$self->{'number_of_atoms_list'} <= 30 || die "Invalid count line" ;
			$self->{'obsolete'} = shift @fields ;
			$self->{'chiral_flag'} = shift @fields ;
			$self->{'number_of_stext_entries'} = shift @fields  ;
			$self->{'number_of_reaction_components'} = shift @fields  ;
			$self->{'number_of_reactants'} = shift @fields  ;
			$self->{'number_of_products'} = shift @fields  ;
			$self->{'number_of_intermediates'} = shift @fields  ;
			$self->{'number_of_add_prop_lines'} = shift @fields  ;
			$self->{'ctab_version'}  = shift @fields ;
		}}
			


		$self->{'original_count_line'} = $_ ; #used by CPSS		
				
			
		#read  blocks for atoms , bonds , ...  
		$self->read_blocks();
		
		# Reading the properties block
		#  lines  like 'M  CHG  6  15   1  16  -1  17   1  18  -1  19   1  20  -1'
		#Two cases:
			# there is a version stamp (v2000 , v3000, v...)
				#then read until M END
			# there is no version stamp
				#then read the correct number of entries
		my @properties_array = () ;
		my $number_of_add_prop_lines = 0 ;
		
		if( $self->version_stamp_is_present())
		{
			$number_of_add_prop_lines = 10e38 ; #infinite number
		}
		else
		{			
			if( defined $self->{'number_of_add_prop_lines'} )
			{
				$number_of_add_prop_lines = $self->{'number_of_add_prop_lines'} ;
			}
		}
		
		my $i ;  #Note for(1 .. 10e38) does not work
		for( $i = 0  ; $i <  $number_of_add_prop_lines; $i++)
		{
			$_ = <STDIN> ;
			defined $_ || die "Insufficient additional properties" ;
			chomp ;

			if( /^\$\$\$\$/ )
			{
				die "Error : unexpected \$\$\$\$";
			}			
			push( @properties_array, $_ ) ;
			
			if( /^M  END/) { last ; }

			 # patch 27-4-96 for reading daselect.rdf
			if( /^\$DATM / ) { last ;}			
	
		} #end for

		#Checking:
		if( ! $self->version_stamp_is_present())
		{		
			if( $number_of_add_prop_lines != $#properties_array)
			{
				print STDERR "$number_of_add_prop_lines $#properties_array\n";
				warn "The number of propertie lines in the counts line does not match the number of properties read";
			}
		}
		
		$self->{'properties_block'} = \@properties_array ;
		
		#$self->inspect() ;
		$self; 		

}
#--------------------------------------------------------------------------#
#Read atoms , bonds, atom-list and Stext blocks
sub read_blocks
{
	my $self = shift ;
	my $mol_string = "" ;
	
	my $number_of_lines_to_be_read = $self->{'number_of_atoms'} + $self->{'number_of_bonds'} 
		+ $self->{'number_of_atoms_list'} + $self->{'number_of_stext_entries'} ;

	for( 1  .. $number_of_lines_to_be_read) {
	    $_ =  <STDIN> ; defined $_ || die "File format error: not enough data";
		$mol_string .= $_ ;
	}
	
	$self->{'mol_string'} = $mol_string ;
		
	$self ;
}

#--------------------------------------------------------------------------#
sub write	{
	my $self = shift ;

	local $_ ;
	
		foreach $line (qw(line1 line2 line3))
		{
			print($self->{$line}, "\n") || die "Write error" ;
		}
	
		$_ = sprintf( "%3d%3s%3s%3s%3s%3s%3s%3s%3s%3s%3s%6s\n",
			$self->{ 'number_of_atoms'},
			$self->{ 'number_of_bonds'},
			$self->{'number_of_atoms_list'},
			$self->{'obsolete'},
			$self->{'chiral_flag'},
			$self->{'number_of_stext_entries'},
			$self->{'number_of_reaction_components'},
			$self->{'number_of_reactants'},
			$self->{'number_of_products'},
			$self->{'number_of_intermediates'},
			$self->{'number_of_add_prop_lines'},
			$self->{'ctab_version'}
		);
		
		#Remove ending blanks 
		s/\s+$// ;
		print( $_, "\n" )|| die "Write error";
			
			
	$self->write_blocks() ;
	

	foreach $prop( @{$self->{'properties_block'}} )	{
		print( "$prop\n") || die "Write error";
	}
	
	$self ;

}

#--------------------------------------------------------------------------#
#Read atoms , bonds, atom-list and Stext blocks, which were saved in the mol_string variable
#See read_blocks
sub write_blocks
{
	my $self = shift ;
	
	defined($self->{'mol_string'} ) || die "Assertion failed" ;
	print( $self->{'mol_string'}) || die "Write error";
	
	$self ;

}

#**************************************************************************#
package MDL_molecule ;
#**************************************************************************#

	@ISA = qw( MDL_non_parsed_molecule );

	@fields_names = (
	 )
	  ;

	
	eval( ${HashObject'new_SubClass_String} ) ;

#--------------------------------------------------------------------------#
#Read atoms , bonds, atom-list and Stext blocks
sub read_blocks
{
	my $self = shift ;

	my @atoms = () ;
	for( $j = 1 ; $j <= $self->{'number_of_atoms'} ; $j ++ ) {
	   	$_ =  <STDIN> ; defined $_ || die "File format error: not enough data";
		push( @atoms, readFromCurrentLine MDL_atom) ;
	}
	$self->{'atoms'} = \@atoms ; #store the references

	my @bonds = () ;
	for( 1  .. $self->{'number_of_bonds'}) {
	    $_ =  <STDIN> ; defined $_ || die "File format error: not enough data";
		push( @bonds, readFromCurrentLine MDL_bond) ;
	}
	
	$self->{'bonds'} = \@bonds ;

	my @list = () ;
	
	for( 1 .. $self->{'number_of_atoms_list'})
	{
	    $_ =  <STDIN> ; defined $_ || die "File format error: not enough data";
		push( @list, readFromCurrentLine MDL_atom_list_line) ;
	}
	$self->{'atom_list_block'} = \@list ;
	
	my @list2 = () ;
	
	for( 1 .. $self->{'number_of_stext_entries'})
	{
	    $_ =  <STDIN> ; defined $_ || die "File format error: not enough data";
		chomp ;
		push( @list2, $_) ;
	}
	$self->{'stext_block'} = \@list2 ;
		
	$self ;
}

#--------------------------------------------------------------------------#
sub write_blocks
{
	my $self = shift ;
	
	foreach $line ( 
		@{$self->{'atoms'}},
		@{$self->{'bonds'}},
		@{$self->{'atom_list_block'}} )
	{
		$line->writeLine() ;
	}
	
	
	foreach $line ( @{$self->{'stext_block'}})
	{
		print( $line, "\n") || die "Write error" ;
	}
	
	$self ;


}
#--------------------------------------------------------------------------#
#Is this atom connected to O N or S?

sub isThisAtomAPolarHydrogen
{
	my $self = shift;
	my $atom = shift ;
	
	$atom->{'symbol'} eq 'H' || return 0 ;
	
	my @connections = $self->connections_indices($self->index_of_atom($atom));
	$#connections <= 1 || die "Assertion failed: H is connected to more that 1 atom" ;
	$#connections >= 1 || die "Assertion failed: H is not connected" ;
	
	my @atoms = @{$self->{'atoms'}} ;
	defined @atoms || die "assertion failed" ;
	
	my $other_atom = ($atoms[$connections[1]])->{'symbol'};

	#Add or remove  atom type if you wish
	$other_atom =~ m/^[ONSFIP]$/i || 
		$other_atom =~ /^(Cl)|(Br)|(Se)$/i;
	
}
#--------------------------------------------------------------------------#
#Has this atom stereo information?
sub isThisAtomStereo
{
	my $self = shift;
	my $atom = shift ;
		
	my $connections = $self->connections_bonds($self->index_of_atom($atom));

	foreach $bond (@$connections)
	{
		defined $bond || die "Assertion failed" ;
		defined $bond->{'stereo'} || die "assertion failed";
		
		if( $bond->{'stereo'} != 0 )
		{
			return 1 ;
		}
	}
	
	return 0;

}

#--------------------------------------------------------------------------#
#             Check Mapping
#   check_for_duplicate_map_numbers (except 0)
#--------------------------------------------------------------------------#
sub has_duplicated_map_numbers {
	
	my $self = shift ;
	my @check_aa_map = () ;
	
	my $aa_map_number ;
	
	@{$self->{'atoms'}} || die "Undefined atoms" ;
	
	foreach $atom (@{$self->{'atoms'}}) {
		$atom || die "Undefined atom" ;
		$aa_map_number = $atom->{'atom_atom_mapping_number'} ;
		$aa_map_number != 0 || next ;
		$check_aa_map[ $aa_map_number ] ++ ;
		
		if( $check_aa_map[ $aa_map_number ] > 1) {
			return $aa_map_number ;
		} 	
	}
	
	0 ;
}
#--------------------------------------------------------------------------#

sub add_atom	{

		my $self = shift ;
		my $atom = shift ;
		my @atoms = () ;
		
		if(! defined $self->{'atoms'})	{
			$self->{'atoms'} = \@atoms ; #store the references
		} 		
		push( @{$self->{'atoms'}}, $atom ) ;
		$self->{ 'number_of_atoms'} ++ ;
		
		$atom ;
}
#--------------------------------------------------------------------------#

sub remove_atom	{
	my $self = shift ;
	my $atom_to_remove = shift ;
	
	return $self->remove_atom_index($self->index_of_atom($atom_to_remove)) ;

}
#--------------------------------------------------------------------------#

sub remove_atom_index	{

		my $self = shift ;
		my $remove_index = shift ;
		my @atoms = () ;
		my @new_atoms = () ;
		my $found = 0 ;
		
		($remove_index > 0 && $remove_index <= $self->{'number_of_atoms'}) || die "Assertion failed" ;
		@atoms = @{$self->{'atoms'}};
		defined $atoms[$remove_index] || die "Assertion failed" ;
		
		if(! defined $self->{'atoms'})	{
			die "Assertion failed: \@atoms is undefined" ;
		}
		@atoms =  @{$self->{'atoms'}} ;
		foreach $each_atom (  @atoms )
		{
			if( $each_atom == $atom)
			{
				$found ++ ;
				next ;
			}
			push( @new_atoms, $each_atom ) ;
		}
		
		$found == 1 || die "Assertion failed";			

		
#2nd special hydrogen treatment: valence
#
	defined $atoms[$remove_index] || die "Assertion failed" ;
#	
#	if( ($atoms[$remove_index])->{'symbol'} eq 'H')
#	{
#		my @connections = $self->connections_indices($remove_index);
#		$#connections == 1 || die "Assertion failed";
#		my $atom =  $atoms[$connections[1]];
#		defined $atom || die "Assertion failed";
#		$atom->{'valence'} = $atom->{'valence'} + 1 ; #++ does not work
#		
#	}

# 3rd the bonds
	my @bonds  = @{$self->{'bonds'}} ;
	my @new_bonds = () ;
	
	
	foreach $bond (@bonds) {
		$bond->{'first_atom'} == $remove_index && next ;
		$bond->{'second_atom'} == $remove_index && next ;
				
		if($bond->{'first_atom'} > $remove_index ) {
			$bond->{'first_atom'} -- ;		
		}
		if($bond->{'second_atom'} > $remove_index ) {
			$bond->{'second_atom'} -- ;		
		}
		#$bond->inspect() ;
		push( @new_bonds, $bond) ; 
	}
	
	
	#update structure
	$self->{'atoms'} = \@new_atoms ;
	$self->{ 'number_of_atoms'} = $#new_atoms ;
	$self->{'bonds'} = \@new_bonds ;
	$self->{'number_of_bonds'} = $#new_bonds ;
	
	$self;
}
#--------------------------------------------------------------------------#
sub add_bond	{

		my $self = shift ;
		my $bond = shift ;
		my @bonds = () ;
		
		if(! defined $self->{'bonds'})	{
			$self->{'bonds'} = \@bonds ; #store the references
		}
		
		push( @{$self->{'bonds'}}, $bond ) ;
		$self->{ 'number_of_bonds'} ++ ;
		
		$bond ;
}
#--------------------------------------------------------------------------#
#  index_of_atom

sub index_of_atom
{
	my $self = shift ;
	my $atom = shift ;
	
	my $index = 0 ; 
	
	foreach $each (@{$self->{'atoms'}})
	{
		$index ++ ;
		if( $each == $atom )
		{
			return $index ;
		}
	}
	
	return -1 ;#impossible value
}
#--------------------------------------------------------------------------#
#answer a collection containing all the other atom indices with which this atom
#index is connected
sub connections_indices
{
	my $self = shift ;
	my $this_index = shift ;
	my @indices = () ;
	
	($this_index >= 1 && $this_index <= $self->{'number_of_atoms'} ) || die "Assertion failed" ;
	
	foreach $bond (@{$self->{'bonds'}})
	{
		if($bond->{'first_atom'} == $this_index )
		{
			push(@indices, $bond->{'second_atom'}) ;
			next ;
		}
		
		if($bond->{'second_atom'} == $this_index )
		{
			push(@indices, $bond->{'first_atom'});
			next ;
		}
		
	}
	
	return @indices ;	
}
#--------------------------------------------------------------------------#
#answer a collection containing all the  bonds which involve this atom
sub connections_bonds
{
	my $self = shift ;
	my $this_index = shift ;


	my @connecting_bonds = () ;
	
	($this_index >= 1 && $this_index <= $self->{'number_of_atoms'} ) || die "Assertion failed" ;
	
	foreach $bond (@{$self->{'bonds'}})
	{
		if($bond->{'first_atom'} == $this_index || $bond->{'second_atom'} == $this_index)
		{
			push(@connecting_bonds, $bond) ;
		}
		
	}
	
	return \@connecting_bonds ;
}
#--------------------------------------------------------------------------#

#       copy_without_atom

# answer a copy of self with one atom removed
# arg : one atom
#--------------------------------------------------------------------------#

sub copy_without_atom
{
	my $self = shift ;
	my $atom = shift ;
	
	defined $atom || die "Assertion failed" ;
	
	my $index = $self->index_of_atom($atom) ;
	$index > 0 || die "Assertion failed" ;
	return $self->copy_without_atom_index($index);
}




#--------------------------------------------------------------------------#

#       copy_without_atom_index

# answer a copy of self with one atom removed
# arg : the index of the atom
#--------------------------------------------------------------------------#

sub copy_without_atom_index
#not tested
{
	my $self = shift ;
	my $remove_index = shift ;
	
	my $copy = $self->copy ;
	
	$remove_index > 0 && $remove_index <= $self->{'number_of_atoms'} || die "Assertion failed" ;


# 1st , the atoms
	my @atoms = @{$self->{'atoms'}} ;
	my @new_atoms = () ;

	$self->{'number_of_atoms'} == $#atoms || die "Assertion failed" ;

	$copy->remove_atom_index( $remove_index );
# 3rd the propertie block
	@{$copy->{'properties_block'}} = @{$self->{'properties_block'}} ; #deep copy
# TO BE COMPLETED	
	



# end : return the copy
	$copy ;	
}

#--------------------------------------------------------------------------#
#Calc geometrical center
#Answer an array ref with the center
sub geometric_center
{
	my $self = shift ;

	my @atoms = @{$self->{'atoms'}} ;

	$self->{'number_of_atoms'} == $#atoms || die "Assertion failed" ;

	my ($x, $y, $z) = (0, 0, 0);
	
	foreach $atom (@atoms)
	{		
		defined( $atom ) || die "assertion failed" ;
		
		$x += $atom->{'x'};	
		$y += $atom->{'y'};
		$z += $atom->{'z'};
		
	}
	
	$x /= $#atoms ;
	$y /= $#atoms ;
	$z /= $#atoms ;
	
	my @a = ($x, $y, $z) ;
	\@a;

	
}
#--------------------------------------------------------------------------#
#Center the molecule using the geometrical center
sub center_coordinates
{
	my $self = shift ;

	my @atoms = @{$self->{'atoms'}} ;

	$self->{'number_of_atoms'} == $#atoms || die "Assertion failed" ;

	my ($cx, $cy, $cz) = @{$self->geometric_center()};
	
	foreach $atom (@atoms)
	{		
		defined( $atom ) || die "assertion failed" ;
		
		$atom->{'x'} -= $cx;	
		$atom->{'y'} -= $cy;	
		$atom->{'z'} -= $cz;	
		
	}
	
	$self ;
	
}
#--------------------------------------------------------------------------#
#translate the molecule 
sub translate_coordinates
{
	my $self = shift ;
	my ($cx, $cy, $cz) = @_;

	defined $cx && defined $cy && defined $cz || die "Assertion failed";
	
	my @atoms = @{$self->{'atoms'}} ;

	$self->{'number_of_atoms'} == $#atoms || die "Assertion failed" ;

	
	foreach $atom (@atoms)
	{		
		defined( $atom ) || die "assertion failed" ;
		
		$atom->{'x'} += $cx;	
		$atom->{'y'} += $cy;	
		$atom->{'z'} += $cz;	
		
	}
	
	$self ;
	
}
#--------------------------------------------------------------------------#
#Scale the molecule using a 3D array
sub scale_coordinates
{
	my $self = shift ;
	my ($sx, $sy, $sz) = @_ ;
	
	defined $sx && defined $sy && defined $sz || die "Assertion failed";

	my @atoms = @{$self->{'atoms'}} ;

	$self->{'number_of_atoms'} == $#atoms || die "Assertion failed" ;

	
	foreach $atom (@atoms)
	{		
		defined( $atom ) || die "assertion failed" ;
		
		$atom->{'x'} *= $sx;	
		$atom->{'y'} *= $sy;	
		$atom->{'z'} *= $sz;	
		
	}
	
	$self ;
	
}
#--------------------------------------------------------------------------#
#Zeroing the coordinates 
sub zeroing_coordinates
{
	my $self = shift ;
	

	my @atoms = @{$self->{'atoms'}} ;

	$self->{'number_of_atoms'} == $#atoms || die "Assertion failed" ;

	
	foreach $atom (@atoms)
	{		
		defined( $atom ) || die "assertion failed" ;
		
		$atom->{'x'} = 0;	
		$atom->{'y'} = 0;	
		$atom->{'z'} = 0;	
		
	}
	
	$self ;
	
}	
#--------------------------------------------------------------------------#
sub remove_stereo_information
{
	my $self = shift ;
	
	$self->zeroing_atom_stereo_parity() ;
	$self->remove_stereo_information_in_bonds() ;

	
	$self ;
	
}	
#--------------------------------------------------------------------------#
#Remove  atom_stereo_parity
sub zeroing_atom_stereo_parity
{
	my $self = shift ;
	

	my @atoms = @{$self->{'atoms'}} ;

	$self->{'number_of_atoms'} == $#atoms || die "Assertion failed" ;

	
	foreach $atom (@atoms)
	{		
		defined( $atom ) || die "assertion failed" ;
		
		$atom->{'atom_stereo_parity'} = 0;			
	}
	
	$self ;
	
}#--------------------------------------------------------------------------#

sub remove_stereo_information_in_bonds
{
	my $self = shift ;
		
	foreach $bond	( @{$self->{'bonds'}} )	{
		$bond->{'stereo'} = 0 ;
	}
}
#--------------------------------------------------------------------------#
#Remove  charges
sub zeroing_atom_charges
{
	my $self = shift ;
	

	my @atoms = @{$self->{'atoms'}} ;

	$self->{'number_of_atoms'} == $#atoms || die "Assertion failed" ;

	
	foreach $atom (@atoms)
	{		
		defined( $atom ) || die "assertion failed" ;
		
		$atom->{'charge'} = 0;			
	}
	
	$self ;
	
}#--------------------------------------------------------------------------#

#       remove_non_polar_hydrogens
# i.e. connected to a carbon
#.must not be used to define a stereo center
# remove non polar hydrogens
#--------------------------------------------------------------------------#

sub remove_non_polar_hydrogens
{
	my $self = shift ;

	my @atoms = @{$self->{'atoms'}} ;

	$self->{'number_of_atoms'} == $#atoms || die "Assertion failed" ;

	foreach $atom (@atoms)
	{		
		defined( $atom ) || die "assertion failed" ;
		
		if( $atom->{'symbol'} eq 'H' && (! $self->isThisAtomAPolarHydrogen($atom ) && (! $self->isThisAtomStereo($atom))))
		{
			$self->remove_atom($atom) ;
			#$self->inspect();
		}
		
	}
	
	
	$self;
	
}
#--------------------------------------------------------------------------#

#       remove all hydrogens
#--------------------------------------------------------------------------#

sub remove_hydrogens
{
	my $self = shift ;

	my @atoms = @{$self->{'atoms'}} ;

	$self->{'number_of_atoms'} == $#atoms || die "Assertion failed" ;

	foreach $atom (@atoms)
	{		
		defined( $atom ) || die "assertion failed" ;
		
		if( $atom->{'symbol'} eq 'H' )
		{
			$self->remove_atom($atom) ;
			#$self->inspect();
		}
		
	}
	
	$self;
	
}
#-----------------------------------------------------------------------#
sub calc_smiles

# Return a string or undef if an error occured
# Uses babel
{
	my $self = shift ;
	my $opt_correct_nitro = shift ;
	
	$ofile = "$$.mol" ;    #generate Babels 'input  file name
	$ifile = "$$.smiles" ; #generate Babels ' output file name	
	
	$self->writeToOutputFile($ofile) ;

	my $smiles_string;
	
	#See the Perl blue book page 230 for handling the return code
	#NOTE : This is buggy !!!! Ctrl C does not work correctly
	
	my $rc = 0xffff & system "babel -imdl $ofile -osmiles $ifile" ;
	
	if( $rc == 0) #Everything went fine
	{

		open( FILEIN, "$ifile") || die "Cannot open $ifile for reading" ;
		
		local $_ = <FILEIN> ;
		chomp;
		
		if( defined $opt_correct_nitro && $opt_correct_nitro == 1)
		{
			#The most specific substitution must come first
			s/N\(=O\)\(O\)/[N+](=O)([O-])/g;
			s/N\(=O\)O/[N+](=O)[O-]/g;
		}
		
		close FILEIN ;
		
		$smiles_string	= $_ ;
					
	} elsif( $rc == 0xffff)
	{
		warn  " Babel command failed \n" ;
		
	}elsif( $rc> 0x80)
	{
		$rc >>= 8 ;
		
		warn " ran with non zero status $rc \n" ;
		
	}else
	{
		if( $rc & 0x80 )
		{
			$rc = &= ^0x80 ;
			print STDERR "coredumped ! " ;
		}
		
		print STDERR " Signal = $rc \n" ;
	}
	
	
	# remove temporary files
	unlink $ofile ;
	unlink $ifile ;
	
	return $smiles_string;
	

}


#-----------------------------------------------------------------------#
#--------------------------------------------------------------------------#

#       randomize the atoms order in the connection table
#--------------------------------------------------------------------------#

sub randomize_atoms_order
{
	my $self = shift ;


	my @properties = @{$self->{'properties_block'}} ;
	
	foreach $prop_line (@properties)
	{
		local $_ = $prop_line;
		 /^M\s+END/ && next ;
		if( /^M\s+(CHG)|(RAD)|(ISO)/ )
		{
			die "Modifications of the property block has yet to be implemented.";
		}
	
	}
	
	
	my @atoms = @{$self->{'atoms'}} ;

	$self->{'number_of_atoms'} == $#atoms || die "Assertion failed" ;
	
	my $random_array_ref = main'random_reordered_progression($self->{'number_of_atoms'});
	#Note $random_array_ref is calculated with a different  array indexing
	# $[
	
	my @new_atoms = ();
	my $j = 1 ;
	foreach $i (@$random_array_ref)
	{
		my $atom = $atoms[$j] ;
		defined( $atom ) || die "assertion failed" ;
		$new_atoms[ $i + $[  ] = $atom ;
		$j ++ ;
		
	}

	$self->{'atoms'} = \@new_atoms ;

	foreach $bond (@{$self->{'bonds'}})
	{
		defined $bond || die ;
		my $first = $bond->{'first_atom'};
		my $second = $bond->{'second_atom'};
		(defined $first && defined $second) || die"Assertion failed" ;
		 $first !=  $second || die"Assertion failed" ;

		($first > 0 && $first <= $#atoms) || die "Assertion failed" ;
		($second > 0 && $second <= $#atoms) || die "Assertion failed" ;
		
		my $new_first_index = $first  ;
		my $new_first = $random_array_ref->[$new_first_index] + $[;
		
		my $new_second_index = $second  ;
		my $new_second = $random_array_ref->[$new_second_index ] + $[;
#		printf STDERR "%d %d %d %d\n", $[, $first, $second, $#atoms;
#		printf STDERR "%d %d %d %d\n", $[, $new_first, $new_second, $#atoms;
#		printf STDERR "%d %d %d %d\n", $[, $random_array_ref->[$first - $[], $random_array_ref->[$second- $[], $#atoms;

		$bond->{'first_atom'} = $new_first;
		$bond->{'second_atom'} = $new_second ;
		$bond->{'first_atom'} != $bond->{'second_atom'} || die "Assertion failed" ;
		
	}

	$self;
	
}
#--------------------------------------------------------------------------#

sub randomize_bonds_order
{
	my $self = shift ;

	my @bonds = @{$self->{'bonds'}} ;

	$self->{'number_of_bonds'} == $#bonds || die "Assertion failed" ;
	
	my $random_array_ref = main'random_reordered_progression($self->{'number_of_bonds'});
	#Note $random_array_ref is calculated with a different  array indexing
	# $[
	
	my @new_bonds = ();
	my $j = 1 ;
	foreach $i (@$random_array_ref)
	{
		my $bond = $bonds[$j] ;
		defined( $bond ) || die "assertion failed" ;
		$new_bonds[ $i +$[] = $bond ;
		$j ++ ;
		
	}
	$self->{'bonds'} = \@new_bonds ;
}


#**************************************************************************#
package MDL_reaction_components ;
#**************************************************************************#

	@ISA = qw( HashObject );

	@fields_names = ('header_block',  'number_of_reactants', 'number_of_products', 'reactants', 'products') ;

	eval( ${HashObject'new_SubClass_String} ) ;

#--------------------------------------------------------------------------#
	
	sub readFromInput {
		my $class = shift ;
		my $self = new($class) ;
		my $n_reactants ;
		my $n_products ;
		my @reactants = () ;
		my @products = () ;
		

		while( <STDIN> ) {
		#       ISIS     101219941611
			if( !/^\s+$/ ) {
				$self->{'header_block'} = $_ ;
				last ;
			}
		}
		while( <STDIN> ) {
		#  2  1
			if( /^([0-9 ][0-9 ][0-9])([0-9 ][0-9 ][0-9])/ ) {
				$n_reactants = $1 ;
				$n_products = $2 ;
				last ;
			}
		
		}
		
		if( $n_reactants <= 0 || $n_products <= 0 ) {
			die( "Number reactants / products not defined") ;
		}
		
		for( $j = 1 ; $j <= $n_reactants ; $j ++ ) {

			# Change BB 28-3-95	
			while( <STDIN> )	{
				if( /^\$MOL/ )	{ last ;}
			}

			push( @reactants, readFromInput MDL_molecule) ;
		}
			
		for( $j = 1 ; $j <= $n_products ; $j ++ ) {
			# Change BB 28-3-95	
			while( <STDIN> )	{
				if( /^\$MOL/ )	{ last ; }
			}
			push( @products, readFromInput MDL_molecule) ;
		}
	
		
		

		
		$self->{'reactants'} = \@reactants ;
		$self->{'products'} = \@products ;
		$self->{'number_of_reactants'} = $n_reactants ;
		$self->{'number_of_products'} = $n_products ;
	
		
		$self; 		

	}
#--------------------------------------------------------------------------#
 # check mapping 
	sub has_a_molecule_with_duplicated_map_numbers
	{
		my $self = shift ;
		my @all_molecules = (@{$self->{'reactants'}}, @{$self->{'products'}}) ;
		
		foreach $molecule ( @all_molecules)	{
			if( $molecule->has_duplicated_map_numbers() ) {
				return 1 ;
			}
		}
		
		return 0 ;
	}

#--------------------------------------------------------------------------#
 # check mapping 
	sub has_a_reactant_with_duplicated_map_numbers
	{
		my $self = shift ;
		my @molecules = (@{$self->{'reactants'}}) ;
		
		foreach $molecule ( @molecules)	{
			if( $molecule->has_duplicated_map_numbers() ) {
				return 1 ;
			}
		}
		
		return 0 ;
	}
#--------------------------------------------------------------------------#
 # check mapping 
	sub has_a_product_with_duplicated_map_numbers
	{
		my $self = shift ;
		my @molecules = (@{$self->{'products'}}) ;
		
		foreach $molecule ( @molecules)	{
			if( $molecule->has_duplicated_map_numbers() ) {
				return 1 ;
			}
		}
		
		return 0 ;
	}


#--------------------------------------------------------------------------#
	# It seems that rulemaker1 does not like stereo information !
	
	sub remove_stereo_information_in_bonds	{
		my $self = shift ;
		my @all_molecules = (@{$self->{'reactants'}}, @{$self->{'products'}}) ;
		
		foreach $molecule ( @all_molecules)	{
			$molecule->remove_stereo_information_in_bonds() ;
		}
	}

	
1 ;	

