view release on metacpan or  search on metacpan

examples/g09/g09_out_pdb.pl  view on Meta::CPAN

#!/usr/bin/env perl
# DMR April 30, 2014
#
#   perl examples/g09_pdb.pl ~/some/path
#
# pull coordinates (all) and charges from Gaussian output (path submitted
# on commandline)
# write out pdbs in tmp directory with charges in the bfactor column..
#

use Modern::Perl;
use HackaMol;
use HackaMol::X::Calculator;
use Math::Vector::Real;
use Path::Tiny;
use File::chdir;

my $path = shift || die "pass path to gaussian outputs";

my $hack = HackaMol->new( data => $path, );

foreach my $out ( $hack->data->children(qr/\.out$/) ) {

    my $Calc = HackaMol::X::Calculator->new(
        mol     => HackaMol::Molecule->new,
        out_fn  => $out,
        map_out => \&output_map,
        scratch => 'tmp',
    );

    local $CWD = $Calc->scratch;
    my $pdb = $Calc->out_fn->basename;
    $pdb =~ s/\.out/\.pdb/;

    $Calc->map_output;
    my $mol = $Calc->mol;

    $_->bfact( $_->charge ) foreach $mol->all_atoms;
    $mol->print_pdb_ts([0 .. $mol->tmax], $pdb);

}

#  our function to map molec info from output
sub output_map {
    my $calc  = shift;
    my $resn  = shift || "TMP";
    my $resid = shift || 1;
    my @lines = $calc->out_fn->lines;

    #my @qs    = nbo_qs(@lines);
    my @qs    = mulliken_qs(@lines);
    my @atoms = Zxyz(@lines);

    die "number of charges not equal to number of atoms" if ( @qs != @atoms );

    #add info for pdb printing
    my $i = 1;
    foreach my $at (@atoms) {
        $at->serial($i);
        $at->resname($resn);
        $at->resid($resid);
        $i++;
    }

    $atoms[$_]->push_charges( $qs[$_] ) foreach 0 .. $#qs;
    $calc->mol->push_atoms(@atoms);

}

sub Zxyz {