#*******************************************************************************
# Copyright (C) 2001-2004 by Novartis Institutes of Biomedical Research
# 
# This program is free software; you can redistribute it and/or modify it
# under the terms of the Novartis Open Source License as published by
# Novartis AG; either version 1, or (at your option) any later version of
# the license.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# Novartis Open Source License for more details.
# 
#*******************************************************************************

#!/usr/bin/env perl
use strict;
use POSIX;
use Getopt::Long;
use Pod::Usage;
use Term::ANSIColor;
use Term::ReadLine;
use IO::File;
use File::Temp "tempfile";
use File::Basename;

my $PRG_PATH=`pwd`;
chomp($PRG_PATH);

my $BASEPATH = getProgPath();

my $HOME  = $ENV{HOME};
my $BIN   = "$BASEPATH/../src";
my $SHARE = "$BASEPATH";

my $MSMSBIN     = (exists($ENV{MSMS})) ? $ENV{MSMS} : "msms";
my $SURFCONVERT = (exists($ENV{SURFCONVERT})) ? $ENV{SURFCONVERT} : "surfconvert";
my $PROPGEN     = (exists($ENV{PROPGENERATOR})) ? $ENV{PROPGENERATOR} : "propgenerator";

my %surfsetup = (
								 pdbfile => "",
								 probe => 1.4,
								 density => 3.0,
								 path  => $PRG_PATH,
								 format => "msms",
								 cutoff => 1.0,
								 occupancy => "esp"
								);



# READLINE
my $term = new Term::ReadLine 'Surface Comp';
my $OUT = $term->OUT;
$term->ornaments(0);

#********************************************************************************
# MAIN
#********************************************************************************

#-------------------------------------------------------------------------------
# OPTIONS and checking


my $help;
my $man;

GetOptions('help' => \$help,
					 'man|m' => \$man,
					 'pdb|i=s' => \$surfsetup{pdbfile},
					 'path|o=s' => \$surfsetup{path},
           'format|f=s' => \$surfsetup{format},
					 'probe=f' => \$surfsetup{probe},
					 'density|d=f' => \$surfsetup{density},
					 'cutoff|c=f' => \$surfsetup{cutoff},
					 'occupancy=s' => \$surfsetup{occupancy});

GetOptions() or pod2usage(2);

pod2usage(1) if $help;
pod2usage(-exitstatus => 0, -verbose => 2) if $man;

if (! -f $surfsetup{pdbfile} )
{
	print(STDERR "The pdb-file $surfsetup{pdbfile} does not exist!\n");
	pod2usage(3);
}

if (! -d $surfsetup{path} )
{
	print(STDERR "The given ouput directory $surfsetup{path} does not exist!\n");
	pod2usage(4);
}

my %dirs = checkPath(%surfsetup)
 or
	die ("The experiment data path ".
			 "$surfsetup{path} does not have the sufficient directory structure!");


#-------------------------------------------------------------------------------
# SURFACE AND PROPERTY GENERATION

my @atoms = extractAtomSpheres(pdbfile => $surfsetup{pdbfile});

my $xyzrfile = makeXYZR(atoms => [ @atoms ]);

my $basename = makeSurface( xyzr => $xyzrfile, dest => $dirs{surf}, %surfsetup);
makeProperties( name => $basename, %surfsetup, %dirs );

system("cp $surfsetup{pdbfile} $dirs{pdb}/$basename.pdb");
unlink($xyzrfile);

1;


#*******************************************************************************
sub getProgPath
#*******************************************************************************
{
	return dirname($0);
}

#*******************************************************************************
sub checkPath
#*******************************************************************************
{
	my %params = @_;

	my %dirnames = ( surf   => $params{path}."/".$params{format},
									 data   => $params{path}."/data",
									 pdb    => $params{path}."/pdb",
									 center => $params{path}."/center");

	foreach my $dir (values(%dirnames))
	{
		return 0 if (! -d $dir);
	}

	return %dirnames;
}

#*******************************************************************************
sub extractAtomSpheres
#*******************************************************************************
{
	my %params = @_;
	my $filename = $params{pdbfile};

	my @atoms;

	open(PDBFILE, "<$filename");

	while (<PDBFILE>)
	{
		chomp;
		
		if (/^ATOM/ || /^HETATM/)
		{
			my $data = substr($_, 30);
			$data =~ s/^\s+//;
			my ($x, $y, $z, $occupancy, $vdwradius, @rest) = split(/\s+/, $data);

			push(@atoms, { x => $x, 'y' => $y, z => $z,
										occupancy => $occupancy,
										vdwradius => $vdwradius });
		}
	}

	close(PDBFILE);

	return @atoms;
}

#*******************************************************************************
sub makeXYZR
#*******************************************************************************
{
	my %params = @_;
	my $atoms = $params{atoms};


	open(XYZRFILE, ">temp.xyzr");

	foreach my $atom (@{$atoms})
	{
		printf(XYZRFILE "%12.6f %12.6f %12.6f %12.6f\n", 
					 $atom->{x}, $atom->{'y'}, $atom->{z}, $atom->{vdwradius});
	}

	close(PDBFILE);

	return "temp.xyzr";
}

#*******************************************************************************
sub makeSurface
#*******************************************************************************
{
	my %params = @_;

	my $xyzr    = $params{xyzr};
	my $density = $params{density};
	my $probe   = $params{probe};
	my $format  = $params{format};
	my $dest    = $params{dest};

	my ($base, $path, $suffix) = fileparse($params{pdbfile}, qr{\.pdb});
	my $basename = $base."_MSMS_".sprintf("%.1f",$density);

	system("$MSMSBIN -if $xyzr -of $basename -d $density -p $probe");
	(-f "$basename.vert")
		or die("Could not generate the MSMS surface for $xyzr");
	
	system("cat $basename.vert $basename.face > $basename.msms");
	unlink($basename.'.vert');
	unlink($basename.'.face');

	system("$SURFCONVERT -d -i $basename.msms -o $dest/$basename.$format msms $format");
	unlink($basename.'.msms');

	return $basename;
}

#*******************************************************************************
sub makeProperties
#*******************************************************************************
{
	my %params = @_;

	my $basename = $params{name};
	my $datadir  = $params{data};
	my $surfdir  = $params{surf};
	my $occprop  = $params{occupancy};
	my $cutoff = $params{cutoff};
	my $format  = $params{format};
	my $pdb     = $params{pdbfile};

	system("$PROPGEN -d -i $surfdir/$basename.$format -f $format".
				 " -o $datadir/$basename.dat".
				 " -DCUTOFF=$cutoff -DPDB=$pdb curvature $occprop");
}

__END__

__END__

=head1 NAME

msmsprepare

=head1 SYNOPSIS

A script that prepares MSMS surfaces and surface data files from PDB files

B<msmsprepare [ --help --man|m ]>
                  B<--pdb|i> I<PDBFILE>
                  B<--path|o> I<DATAPATH>
                  B<--format|f> I<OUTPUT_SURFACE_FILE_FORMAT>
                  B<--probe> I<PROBE_RADIUS> B<--density|d> I<POINT_DENSITY>
                  B<--cutoff|c> I<GLOBAL_CURVATURE_CUTOFF>
                  B<--occupancy> I<OCCUPANCY_PROPERTY>

=head1 OPTIONS

=over 4

=item --help

print out this help page.

=item --man

print out a detailed help page as manpage

=item --path [default=I<.>]

The path where the generated data files will be stored. (see also
L<preparesurfcomp/"INPUT DATA STRUCTURE">)

=item --format [default=I<msms>]

The file format of the generated surface file. Possible values are off, msms
and sld. Please note that the specification of the format determines the name
of the directory where the surface files are written.

=item --probe [default=I<1.4 Angs.>]

The size of the probe sphere that was used to create the molecular
surfaces. This parameter is for information purposes only and can be left as is
or set to any value if no solvent accessible surfaces have been used.

=item --density|d [default=I<3/ square Angs.>]

The point density of the generated surface file

=item --cutoff|c [default=I<1.0>]

The range cutoff applied in the neighborhood selection of the global curvature
calculation. (see also [Zachmann1992])

=item --occupancy [default=I<esp>]

The meaning of the occupancy value in the pdb file. The default value is B<esp>
which means that the occupancy field should contain partial charges for each
atom.

=back

=head1 DESCRIPTION

The B<msmsprepare> script takes a special formed PDB file, generates molecular
surfaces with the program MSMS [Sanner1995] and calculates various property
values for all the points on the surface. The script writes the files into a
data structure which is appropriate for the L<preparesurfcomp(1)> script.

The program expects a special PDB file where a property value and the atoms'
van der Waals radii are stored for each atom directly after the x, y and z
coordinate at the places of the I<occupancy> and I<tempFactor> variables. The
radii are taken as input for the MSMS program while the values stored in the
I<occupancy> column are taken as a distinct atom-spectifc chemical property. If
you use the default (I<esp>) the occupance column (right after the z
coordinate) should contain the atomic point charges.

=head2 PROPERTY-CALCULATION

At the moment the following properties can be calculated for the molecular
surfaces:

=over 4

=item global/local curvature (CC1, CC2, STI)

The global and local curvature is calculated according to the algorithm given
by Zachmann et. al. and uses a B<cutoff> to select only a distinct neighborhood
around a point for the curvature calculation.

=item electrostatic potential (ESP)

Calculates the electrostatic potential using Coloumbs law. To get reasonable
results the atomic point charges should be given as the I<occupancy> values of
the PDB file.

=head1 REFERENCES

M.F. Sanner et. al., I<Fast and Robust Computation of Molecular Surfaces>,
Proc. 11th ACM Symp. Comp. Geom., B<C6-C7>, 1995.

C.D. Zachmann et. al., I<Topological Analysis of Complex Molecular Surfaces>,
J. Comp. Chem., B<13:1>, 76-84, 1992.
