view release on metacpan or  search on metacpan

examples/g09/g09_in_ss-set.pl  view on Meta::CPAN

#!/usr/bin/env perl
# Demian Riccardi, June 3, 2014
#
# This example takes an xyz/pdb file of a molecule with a disulfide 
# (or modified disulfide R-S-Hg-S-R), rotates the R-S...S-R to a  
# set angle passed at command line, and  generates Gaussian 09 inputs 
# for the B3PW91/[SDD/]6-31+G** level of theory.
# 
use Modern::Perl;
use HackaMol;
use HackaMol::X::Calculator;
use lib 'basis_sets';
use YAML::XS;
use Path::Tiny;

###############################################################################
#          load in the molecule and initialize charge and multiplicity        #
###############################################################################
my $bldr = new HackaMol;
my $file = path(shift);                           # or die "pass xyz/pdb file";
my $dang = shift or die "pass target dihedral";
my $name = $file->basename(qr/\..*/);
my $mol  = $bldr->read_file_mol($file);
$mol->multiplicity(1);
$mol->push_charges(0);
$mol->fix_serial(1);
my @atoms = $mol->all_atoms;

##############################################################################
#                 Set up CS -- SC dihedral                                   #
# Designed for molecules with only two S atoms in disulfide or separated by  #
# atom in between:                                                           #
#    R-C-S.[x].S-C-R where [x] is optional atom                              #
##############################################################################

my @Ss = grep { $_->symbol eq 'S' } @atoms;
my @Cs = grep { $_->symbol eq 'C' } @atoms;

#find S-C bonds
my @SCs = $bldr->find_bonds_brute(
    bond_atoms => [@Ss],
    candidates => [@Cs],
    fudge      => 0.45,
);

#bond_atoms (S) are first in the group! wanted: C-S -- S-C
my ($dihe) =
  $bldr->build_dihedrals( reverse( $SCs[0]->all_atoms ), $SCs[1]->all_atoms );

##############################################################################
#          Find atoms to rotate about dihedral for scan: qrotatable          #
##############################################################################

my $init = {
    $dihe->get_atoms(1)->iatom => 1,
    $dihe->get_atoms(2)->iatom => 1,
};
my $atoms_rotate = qrotatable( $mol->atoms, $dihe->get_atoms(2)->iatom, $init );
delete( $atoms_rotate->{ $Ss[0]->iatom } );
delete( $atoms_rotate->{ $Ss[1]->iatom } );

my $group_rotate =
  HackaMol::AtomGroup->new( atoms => [ @atoms[ keys %{$atoms_rotate} ] ] );

##############################################################################
#          Set up basis sets and input parameters for G09 input              #
##############################################################################
my $opt    = 'opt=modredun freq';    # set to space if single point
my $thr    = 'b3pw91';