package Demeter::Atoms; =for Copyright . Copyright (c) 2006-2019 Bruce Ravel (http://bruceravel.github.io/home). All rights reserved. . This file is free software; you can redistribute it and/or modify it under the same terms as Perl itself. See The Perl Artistic 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. =cut use autodie qw(open close); use Moose; extends 'Demeter'; with 'Demeter::Tools'; with 'Demeter::Atoms::Absorption'; with 'Demeter::Atoms::Cif' if $Demeter::STAR_Parser_exists; use Demeter::StrTypes qw( Element Edge AtomsLattice AtomsGas AtomsObsolete SpaceGroup FileName Empty ); use Demeter::NumTypes qw( Natural PosInt PosNum NonNeg OneToFour FeffVersions ); #use diagnostics; use Carp; use Chemistry::Elements qw(get_Z); use File::Basename; use List::Util qw(min max reduce); use List::MoreUtils qw(any); #use Math::Cephes::Fraction qw(fract); use POSIX qw(ceil); use Safe; use Scalar::Util qw(looks_like_number); use Text::Template; use Xray::Absorption; use Xray::Crystal; use Demeter::Constants qw($NUMBER $SEPARATOR $EPSILON4 $FEFFNOTOK); use Const::Fast; const my $FRAC => 100000; const my %EDGE_INDEX => (k =>1, l1=>2, l2=>3, l3=>4, m1=>5, m2=>6, m3=>7, m4=>8, m5=>9, n1=>10, n2=>11, n3=>12, n4=>13, n5=>14, n6=>15, n7=>16, ); #has 'cell' => (is => 'rw', isa =>Empty.'|Xray::Crystal::Cell', default=> q{}); has 'cell' => (is => 'rw', isa =>'Any', default=> sub{Xray::Crystal::Cell->new;}, ); has 'message_buffer' => (is => 'rw', isa => 'Str', default => q{}); has 'space' => (is => 'rw', isa => 'Str', default => sub{q{}}, trigger => sub{ my ($self, $new) = @_; return if not $new; $new = $self->colon_in_group($new); $self -> cell -> space_group($new); $self->is_populated(0); $self->absorption_done(0); $self->mcmaster_done(0); $self->i0_done(0); $self->self_done(0); }); has $_ => (is => 'rw', isa => 'Bool', default => 0,) foreach (qw(is_rhomb is_hex is_first is_second)); has 'a' => (is => 'rw', isa => NonNeg, default=> 0, trigger => sub{ my ($self, $new) = @_; return if not $new; $self->is_populated(0); $self->absorption_done(0); $self->mcmaster_done(0); $self->i0_done(0); $self->self_done(0); }); has 'b' => (is => 'rw', isa => NonNeg, default=> 0, trigger => sub{ my ($self, $new) = @_; return if not $new; $self->is_populated(0); $self->absorption_done(0); $self->mcmaster_done(0); $self->i0_done(0); $self->self_done(0); }); has 'c' => (is => 'rw', isa => NonNeg, default=> 0, trigger => sub{ my ($self, $new) = @_; return if not $new; $self->is_populated(0); $self->absorption_done(0); $self->mcmaster_done(0); $self->i0_done(0); $self->self_done(0); }); has 'alpha' => (is => 'rw', isa => NonNeg, default=> 90, trigger => sub{ my ($self, $new) = @_; return if not $new; $self->is_populated(0); $self->absorption_done(0); $self->mcmaster_done(0); $self->i0_done(0); $self->self_done(0); }); has 'beta' => (is => 'rw', isa => NonNeg, default=> 90, trigger => sub{ my ($self, $new) = @_; return if not $new; $self->is_populated(0); $self->absorption_done(0); $self->mcmaster_done(0); $self->i0_done(0); $self->self_done(0); }); has 'gamma' => (is => 'rw', isa => NonNeg, default=> 90, trigger => sub{ my ($self, $new) = @_; return if not $new; $self->is_populated(0); $self->absorption_done(0); $self->mcmaster_done(0); $self->i0_done(0); $self->self_done(0); }); has 'rmax' => (is => 'rw', isa => NonNeg, default=> sub{ shift->co->default("atoms", "rmax") || 8}, trigger => sub{ my ($self, $new) = @_; $self->is_expanded(0) if $new}); has 'rpath' => (is => 'rw', isa => NonNeg, default=> sub{ shift->co->default("atoms", "rpath") || 5}, trigger => sub{ my ($self, $new) = @_; $self->is_expanded(0) if $new}); has 'rscf' => (is => 'rw', isa => NonNeg, default=> sub{ shift->co->default("atoms", "rscf") || 5},); has 'do_scf' => (is => 'rw', isa =>'Bool', default=> 1); has 'do_xanes' => (is => 'rw', isa =>'Bool', default=> 0); has 'rss' => (is => 'rw', isa => NonNeg, default=> 0); has 'edge' => (is => 'rw', isa => Empty.'|'.Edge, coerce => 1, default=> q{}, trigger => sub{ my ($self, $new) = @_; if (exists($EDGE_INDEX{lc($new)})) { my ($central, $xcenter, $ycenter, $zcenter) = $self -> cell -> central($self->core); $self->iedge($EDGE_INDEX{lc($new)}); $self->eedge(Xray::Absorption->get_energy($central->element, $new)) if ($central =~ m{Site}); } else { $self->iedge(0); $self->eedge(0); }; }); has 'iedge' => (is => 'rw', isa => Natural, default=> 1); has 'eedge' => (is => 'rw', isa => NonNeg, default=> 0); has 'core' => (is => 'rw', isa =>'Str', default=> q{}); has 'corel' => (is => 'rw', isa =>'Str', default=> q{}); has 'partial_occupancy' => (is => 'rw', isa =>'Bool', default=> 0); has 'shiftvec' => ( traits => ['Array'], is => 'rw', isa => 'ArrayRef', default => sub { [0, 0, 0] }, handles => { 'push_shiftvec' => 'push', 'pop_shiftvec' => 'pop', 'clear_shiftvec' => 'clear', } ); has 'polarization' => ( traits => ['Array'], is => 'rw', isa => 'ArrayRef', default => sub { [0, 0, 0] }, handles => { 'push_polarization' => 'push', 'pop_polarization' => 'pop', 'clear_polarization' => 'clear', } ); has 'ellipticity' => ( traits => ['Array'], is => 'rw', isa => 'ArrayRef', default => sub { [0, 0, 0] }, handles => { 'push_ellipticity' => 'push', 'pop_ellipticity' => 'pop', 'clear_ellipticity' => 'clear', } ); has 'file' => (is => 'rw', isa =>FileName, default=> q{}, trigger => sub{ my ($self, $new) = @_; if ($new) { $self->read_inp; #$self->update_edge; }; }); has 'cif' => (is => 'rw', isa =>FileName, default=> q{}, trigger => sub{ my ($self, $new) = @_; if ($new) { if (not $Demeter::STAR_Parser_exists) { warn "STAR::Parser is not available, so CIF files cannot be imported"; return; }; $self->read_cif; #$self->update_edge; }; }); has 'record' => (is => 'rw', isa => NonNeg, default=> 0, trigger => sub{ my ($self, $new) = @_; if (not $Demeter::STAR_Parser_exists) { warn "STAR::Parser is not available, so CIF files cannot be imported"; return; }; $self->read_cif if ($self->cif); }); has 'titles' => ( traits => ['Array'], is => 'rw', isa => 'ArrayRef[Str]', default => sub { [] }, handles => { 'push_titles' => 'push', 'pop_titles' => 'pop', 'clear_titles' => 'clear', } ); has 'ipot_style' => (is => 'rw', isa =>'Str', default=> sub{ shift->mo->config->default("atoms","ipot_style") || 'elements'}, trigger => sub{ my ($self, $new) = @_; $self->is_ipots_set(0) if $new}); has 'feff_version' => (is => 'rw', isa => FeffVersions, default=>sub{ shift->mo->config->default("atoms","feff_version") || 6}); has 'nitrogen' => (is => 'rw', isa => NonNeg, default=> 0, trigger => sub{ my ($self, $new) = @_; ($new) ? $self->gases_set(1) : $self->gases_set(0) }); has 'argon' => (is => 'rw', isa => NonNeg, default=> 0, trigger => sub{ my ($self, $new) = @_; ($new) ? $self->gases_set(1) : $self->gases_set(0) }); has 'xenon' => (is => 'rw', isa => NonNeg, default=> 0, trigger => sub{ my ($self, $new) = @_; ($new) ? $self->gases_set(1) : $self->gases_set(0) }); has 'krypton' => (is => 'rw', isa => NonNeg, default=> 0, trigger => sub{ my ($self, $new) = @_; ($new) ? $self->gases_set(1) : $self->gases_set(0) }); has 'helium' => (is => 'rw', isa => NonNeg, default=> 0, trigger => sub{ my ($self, $new) = @_; ($new) ? $self->gases_set(1) : $self->gases_set(0) }); has 'gases_set' => (is => 'rw', isa =>'Bool', default=> 0); has 'xsec' => (is => 'rw', isa =>'LaxNum', default=> 0); has 'deltamu' => (is => 'rw', isa =>'LaxNum', default=> 0); has 'density' => (is => 'rw', isa =>'LaxNum', default=> 0); has 'mcmaster' => (is => 'rw', isa =>'LaxNum', default=> 0, trigger => sub{ my ($self, $new) = @_; my $n= $self->netsig; $self->netsig($n+$new); }); has 'i0' => (is => 'rw', isa =>'LaxNum', default=> 0, trigger => sub{ my ($self, $new) = @_; my $n= $self->netsig; $self->netsig($n+$new); }); has 'selfamp' => (is => 'rw', isa =>'LaxNum', default=> 0); has 'selfsig' => (is => 'rw', isa =>'LaxNum', default=> 0, trigger => sub{ my ($self, $new) = @_; my $n= $self->netsig; $self->netsig($n+$new); }); has 'netsig' => (is => 'rw', isa =>'LaxNum', default=> 0); has 'is_imported' => (is => 'rw', isa =>'Bool', default=> 0, trigger => sub{ my ($self, $new) = @_; $self->is_populated(0) if ($new==0); }); has 'is_populated' => (is => 'rw', isa =>'Bool', default=> 0, trigger => sub{ my ($self, $new) = @_; $self->is_ipots_set(0) if ($new==0); }); has 'is_ipots_set' => (is => 'rw', isa =>'Bool', default=> 0, trigger => sub{ my ($self, $new) = @_; $self->is_expanded(0) if ($new==0); }); has 'is_expanded' => (is => 'rw', isa =>'Bool', default=> 0); has 'absorption_done' => (is => 'rw', isa =>'Bool', default=> 0); has 'mcmaster_done' => (is => 'rw', isa =>'Bool', default=> 0); has 'i0_done' => (is => 'rw', isa =>'Bool', default=> 0); has 'self_done' => (is => 'rw', isa =>'Bool', default=> 0); has 'sites' => ( traits => ['Array'], is => 'rw', isa => 'ArrayRef', default => sub { [] }, handles => { 'push_sites' => 'push', 'pop_sites' => 'pop', 'clear_sites' => 'clear', } ); has 'cluster' => ( traits => ['Array'], is => 'rw', isa => 'ArrayRef', default => sub { [] }, handles => { 'push_cluster' => 'push', 'pop_cluster' => 'pop', 'clear_cluster' => 'clear', } ); has 'nclus' => (is => 'rw', isa =>'Str', default=> 0); sub BUILD { my ($self, @params) = @_; $self->mo->push_Atoms($self); }; sub DEMOLISH { my ($self) = @_; $self->alldone; }; sub refresh { my ($self) = @_; $self->is_populated(0); }; sub out { my ($self, $key) = @_; my $format = $self->co->default("atoms", "precision") || "9.5f"; $format = '%' . $format; my $val = sprintf("$format", $self->$key); return $val; }; sub clear { my ($self) = @_; $self->$_(0) foreach (qw(a b c nitrogen argon xenon helium krypton gases_set rss eedge)); $self->$_(0) foreach (qw(is_rhomb is_hex is_first is_second do_xanes)); $self->$_(1) foreach (qw(do_scf iedge)); $self->rmax( Demeter->co->default('atoms', 'rmax') ); $self->rpath(Demeter->co->default('atoms', 'rpath')); $self->rscf( Demeter->co->default('atoms', 'rscf') ); $self->$_(90) foreach (qw(alpha beta gamma)); $self->space(q{}); $self->edge(q{}); $self->core(q{}); $self->corel(q{}); $self->clear_sites; $self->clear_cluster; $self->shiftvec([0,0,0]); $self->polarization([0,0,0]); $self->ellipticity([0,0,0]); $self->clear_titles; $self->cell(Xray::Crystal::Cell->new); $self->is_imported(0); $self->is_populated(0); $self->is_ipots_set(0); $self->is_expanded(0); $self->absorption_done(0); $self->mcmaster_done(0); $self->i0_done(0); $self->self_done(0); }; sub read_inp { my ($self) = @_; my $reading_atoms_list = 0; $self->clear; my $file = $self->file; croak("Atoms: no input file provided") if (not $file); croak("Atoms: \"$file\" does not exist") if (not -e $file); croak("Atoms: \"$file\" could not be read") if (not -r $file); #$self->set(file=>$file); open(my $INP, $file); while (my $line = (<$INP>)) { next if ($line =~ m{\A\s*\z}); next if ($line =~ m{\A\s*[\#\%\!\*]}); next if ($line =~ m{\A\s*-{2,}}); chomp $line; $line =~ s{^\s+}{}; $line =~ s{\s+$}{}; ($reading_atoms_list) and do { $self->parse_atoms_line($line); next; }; ($line =~ m{\A\s*title}) and do { $line =~ s{\A\s*title\s*=?\s*}{}; $self->push_titles($line); next; }; ($line =~ m{\A\s*atoms?}) and do { # read the remaining lines as the atoms list $reading_atoms_list = 1; next; }; ## parse each line $self->parse_line($line); }; close $INP; if (not $self->core) { my $maxz = 0; foreach my $site (@{$self->sites}) { my @list = split(/\|/, $site); if (get_Z($list[0]) > $maxz) { $maxz = get_Z($list[0]); $self->core($list[4]); }; }; }; $self->is_imported(1); return $self; }; sub parse_line { my ($self, $line) = @_; #return if not $line; my $file = $self->file; my @words = split(/$SEPARATOR/, $line); my $key = shift @words; (my $rest = $line) =~ s{\A$key$SEPARATOR}{}; if ($key =~ m{space(?:group)?}i) { my $end = (length($rest) < 13) ? length($rest) : 13; my $sg = substr($rest, 0, $end); $self->space($sg); $self->alpha($self->cell->alpha); $self->beta ($self->cell->beta ); $self->gamma($self->cell->gamma); $rest = substr($rest, $end, -1); } else { @words = split(/$SEPARATOR/, $rest); my $val = shift @words; $rest =~ s{$val(?:$SEPARATOR)?}{}; my $vv = ($key =~ m{\bout}) ? shift @words : q{}; $rest =~ s{$vv(?:$SEPARATOR)?}{}; my $vvv = ($key =~ m{shi|daf|qve|ref|pol|ell}) ? shift @words : q{}; $rest =~ s{$vvv(?:$SEPARATOR)?}{}; my $vvvv = ($key =~ m{shi|daf|qve|ref|pol|ell}) ? shift @words : q{}; $rest =~ s{$vvvv(?:$SEPARATOR)?}{}; return if ($key =~ m{\#}); $key = lc($key); my $kk = ($key =~ m{shi}) ? 'shiftvec' : ($key =~ m{pol}) ? 'polarization' : $key; if (($self->meta->has_method($kk)) and ($key =~ m{shi|daf|qve|ref|pol|ell})) { if ((not looks_like_number($val)) or (not looks_like_number($val)) or (not looks_like_number($val))) { carp("\"$key\" takes three numbers as its value\nfound \"$key = $val $vvv $vvvv\"\n($file line $.)\n\n"); return $self; } else { $self->$kk([$val, $vvv, $vvvv]); } } elsif ($self->meta->has_method($key)) { if ((any {$key eq $_} (qw(a b c alpha beta gamma rmax rpath rscf nitrogen argon xenon krypton helium))) and (not looks_like_number($val))) { carp("\"$key\" takes a number as its value\nfound \"$key = $val\"\n($file line $.)\n\n"); return $self; } else { $self->$key(lc($val)); }; } elsif (is_AtomsObsolete($key)) { carp("\"$key\" is a deprecated Atoms keyword ($file line $.)\n\n"); return $self; } else { carp("\"$key\" is not an Atoms keyword ($file line $.)\n\n"); return $self; }; }; $self->parse_line($rest) if (($rest !~ m{\A\s*\z}) and ($rest !~ m{\A\s*[\#\%\!\*]})); return $self; }; sub parse_atoms_line { my ($self, $line) = @_; return 0 if ($line =~ m{\A\s*[\#\%\!\*]}); my ($el, $x, $y, $z, $tag) = split(" ", $line); $tag ||= $el; ($tag = $el) if ($tag =~ m{\A$NUMBER\z}); $tag =~ s{$FEFFNOTOK}{}g; # scrub characters that will confuse Feff my $this = join("|",$el, $x, $y, $z, $tag); $self->push_sites($this); return $self; }; sub colon_in_group { my ($self, $group) = @_; if ($group =~ m{:\s*([12hHrR])\s*\z}) { my $colon = $1; (my $stripped = $group) =~ s{\s*:\s*[12hHrR]\s*\z}{}; $self->is_rhomb(1) if (lc($colon) eq 'r'); $self->is_hex(1) if (lc($colon) eq 'h'); $self->is_first(1) if ($colon eq '1'); $self->is_second(1) if ($colon eq '2'); return $stripped; } else { return $group; }; }; sub populate { my ($self) = @_; my @sites; my $ra = $self->sites; foreach my $s (@$ra) { my ($el, $x, $y, $z, $tag) = split(/\|/, $s); croak("$el is not a valid element symbol\n") if not is_Element($el); next if (lc($el) =~ m{\Anu}); push @sites, Xray::Crystal::Site->new(element=>$el, x=>_interpret($x), y=>_interpret($y), z=>_interpret($z), tag=>$tag); }; ### creating and populating cell return $self if not $self->space; $self -> cell -> space_group($self->space); $self -> handle_colon; foreach my $key (qw(a b c alpha beta gamma)) { my $val = $self->$key; $self -> cell->$key($val) if $val; }; ## Group: $cell->get(qw(given_group space_group class setting)) ## Bravais: $cell->get('bravais') $self -> cell -> shiftvec($self->shiftvec); $self -> cell -> populate(\@sites); foreach my $key (qw(a b c alpha beta gamma)) { $self->$key($self->cell->$key); }; my ($central, $xcenter, $ycenter, $zcenter) = $self -> cell -> central($self->core); $self->update_edge; $self->set(is_populated => 1, corel => ucfirst(lc($central->element)), ); return $self; }; sub element_check { my ($self, $sym) = @_; return is_Element($sym); }; sub _interpret { my ($str) = @_; my $cpt = new Safe; my $retval = $cpt->reval($str); return $retval; }; sub build_cluster { my ($self) = @_; $self->populate if (not $self->is_populated); my ($cell, $core) = $self->get("cell", "core"); my @sites = @{ $cell->sites }; map { $_ -> in_cluster(0) } @sites; #Demeter->pjoin('build_cluster',$self->get(qw(a b c alpha beta gamma))); #Demeter->pjoin('build_cluster',$self->cell->get(qw(txx tyx tyz tzx tzz))); my $rmax = $self->rmax; my @cluster = (); my ($central, $xcenter, $ycenter, $zcenter) = $cell -> central($core); my $setting = $cell->group->setting; my $crystal_class = $cell->group->class; my $do_tetr = ($crystal_class eq "tetragonal" ) && ($setting); #### here my ($aa, $bb, $cc) = $cell -> get("a", "b", "c"); my $xup = ceil($rmax/$aa - 1 + $xcenter); my $xdn = ceil($rmax/$aa - $xcenter); my $yup = ceil($rmax/$bb - 1 + $ycenter); my $ydn = ceil($rmax/$bb - $ycenter); my $zup = ceil($rmax/$cc - 1 + $zcenter); my $zdn = ceil($rmax/$cc - $zcenter); ##print join(" ", "up,dn", $xup, $xdn, $yup, $ydn, $zup, $zdn), $/; #my $num_z = int($rmax/$cc) + 1; # | my $rmax_squared = $rmax**2; # (sprintf "%9.5f", $rmax**2); my ($contents) = $cell -> contents; foreach my $nz (-$zdn .. $zup) { foreach my $ny (-$ydn .. $yup) { foreach my $nx (-$xdn .. $xup) { foreach my $pos (@{$contents}) { my ($x, $y, $z) = ($$pos[1]+$nx, $$pos[2]+$ny, $$pos[3]+$nz); ($x, $y, $z) = ($x-$xcenter, $y-$ycenter, $z-$zcenter); ($x, $y, $z) = $cell -> metric($x, $y, $z); ($do_tetr) and ($x, $y) = (($x+$y)/sqrt(2), ($x-$y)/sqrt(2)); #my ($fx, $fy, $fz) = &rectify_formula(@$pos[4..6], $nx, $ny, $nz); #printf "out: %25s %25s %25s\n\n", $fx, $fy, $fz; my $r_squared = sprintf "%9.5f", $x**2 + $y**2 + $z**2; if ($r_squared < $rmax_squared) { my $this_site = [$x, $y, $z, $$pos[0], $r_squared, # cache the (sprintf "%11.7f", $x), # stuff needed (sprintf "%11.7f", $y), # for sorting (sprintf "%11.7f", $z), #$fx, $fy, $fz, ]; $$pos[0] -> in_cluster(1); push @cluster, $this_site; ## (push @neutral, $this_site); }; }; }; }; }; ## =============================== sort the cluster (& neutral clus.) @cluster = sort { ($a->[4] cmp $b->[4]) # sort by distance squared or ... or ($a->[3] cmp $b->[3]) # by tag alphabetically (using string coercion) or ... or ($a->[7] cmp $b->[7]) # by z value or ... or ($a->[6] cmp $b->[6]) # by y value or ... or ($a->[5] cmp $b->[5]) # by x value ## or ## ($ {$b->[3]}->{Host} <=> $ {$a->[3]}->{Host}); # hosts before dopants } @cluster; if ($#cluster > 499) { carp(sprintf('Your cluster has more than 500 atoms, which is the hard-wired limit for Feff6L. Feff6L was run using only the first 500 atoms. You might want to reduce the value of the cluster size (Rmax = %.2f).\n', $self->rmax)); }; if ($#cluster <= 0) { carp sprintf('You have specified crystal data resulting in 0 scattering atoms. Possible reasons include: * One or more lattice constant is very large * The cluster size (Rmax = %.2f) is too small * You have used cartesian coordinates for sites rather than fractional coordinates ', $self->rmax); return $self; }; $self->set(cluster => \@cluster, nclus => $#cluster+1, rss => sprintf('%'.$self->co->default("atoms", "precision"), $cluster[1]->[4]*$self->co->default("atoms", "smallsphere")), ); # ## final adjustment to the formulas, store the formulas for the # ## central atom ... # $keys -> {cformulas} = # [$$r_cluster[0][8], $$r_cluster[0][9], $$r_cluster[0][10]]; # ## ## ... subtract the central atom coordinates from each site ... # ## foreach my $site (reverse(@$r_cluster)) { # ## (@$site[8..10]) = # ## ($$site[8] . " - Xc", $$site[9] . " - Yc", $$site[10] . " - Zc"); # ## }; # ## ## ... and set the central atom to an empty string # ## ($$r_cluster[0][8], $$r_cluster[0][9], $$r_cluster[0][10]) = ("", "", ""); ## if this is a tetragonal crystal in the C or F setting , rotate ## all the coordinates back to the original setting if ($do_tetr) { my ($a, $b) = $cell->get("a", "b"); $cell->b($a*sqrt(2)); $cell->b($b*sqrt(2)); }; $self->is_expanded(1); return $self; }; sub set_ipots { my ($self) = @_; $self->build_cluster if (not $self->is_expanded); if ($#{$self->cluster} == 0) { carp("You have no scattering atoms in the cluster.\n\n"); return -1; }; my ($cell, $how) = $self->get("cell", "ipot_style"); my @sites = @{ $cell->sites }; my $i = 1; my %seen = (); if ($how =~ m{\Ata}) { foreach my $s (@sites) { if ($s->in_cluster) { my $tag = lc($s->tag); $seen{$tag} = $i++ if (not $seen{$tag}); $s -> ipot($seen{$tag}); }; }; } elsif ($how =~ m{\Ael}) { foreach my $s (@sites) { if ($s->in_cluster) { my $el = lc($s->element); $seen{$el} = $i++ if (not $seen{$el}); $s -> ipot($seen{$el}); }; }; } else { ## sites foreach my $s (@sites) { if ($s->in_cluster) { $s -> ipot($i); ++$i }; }; }; ## get the reduced stoichiometry for feff8's potentials list my @count = (0,0,0,0,0,0,0,0); my $top = -999; foreach my $s (@sites) { my $ipot = $s->ipot; $count[$ipot] += $s->in_cell; $top = max($top, $ipot); }; ## get greatest common divisor (thanks to Demeter::Tools for ## "euclid" (which was swiped from Math::Numbers) and List::Util for ## "reduce") if ($self->co->default("atoms", "gcd")) { my $gcd = reduce { ($self->euclid($a,$b))[0] } @count[1..$top]; foreach my $s (1 .. $top) { $count[$s] /= $gcd; }; }; foreach my $s (@sites) { my $ipot = $s->ipot; $s->stoi($count[$ipot]); }; $self->is_ipots_set(1); if (--$i > 7) { my $ii = $i+1; carp("You have $ii unique potentials, but Feff only allows 7.\n\n"); return -1; }; return 0; }; override 'template' => sub { my ($self, $file, $rhash) = @_; my $cell = $self->cell; my $tmpl = File::Spec->catfile(dirname($INC{"Demeter.pm"}), "Demeter", "templates", "atoms", "$file.tmpl"); croak("Unknown Atoms template file -- type $file: $tmpl") if (not -e $tmpl); my $template = Text::Template->new(TYPE => 'file', SOURCE => $tmpl) or croak("Couldn't construct template: $Text::Template::ERROR"); $rhash ||= {}; my $string = $template->fill_in(HASH => {A => \$self, C => \$cell, %$rhash}, PACKAGE => "Demeter::Templates"); $string ||= q{}; $string =~ s{^\s+}{}; # remove leading white space $string =~ s{\n(?:[ \t]+\n)+}{\n}; # regularize white space between blocks of text $string =~ s{\s+$}{\n}; # remove trailing white space $string =~ s{<<->>\n}{}g; # convert newline token into a real newline $string =~ s{<>}{\n}g; # convert newline token into a real newline $string =~ s{<<( *)>>}{$1}g; #} # convert white space token into real white space return $string; }; sub cluster_list { my ($self, $pattern) = @_; $pattern ||= " %9.5f %9.5f %9.5f %d %-10s %9.5f\n"; $self->set_ipots if (not $self->is_ipots_set); my $string = q{}; my @list = @ {$self->cluster }; my $abs = shift @list; # absorber must be ipot 0 $string .= sprintf($pattern, $abs->[0], $abs->[1], $abs->[2], 0, $abs->[3]->tag, sqrt($abs->[4]) ); my %seen; # index tags by shell foreach my $pos (@list) { if (not defined($seen{$pos->[3]->tag})) { $seen{$pos->[3]->tag} = [1, sqrt($pos->[4])]; }; ++$seen{$pos->[3]->tag}->[0] if (sqrt($pos->[4]) - $seen{$pos->[3]->tag}->[1] > $EPSILON4); # increment index if R has increased my $tag = join(".", $pos->[3]->tag, $seen{$pos->[3]->tag}->[0]); $string .= sprintf($pattern, $pos->[0], $pos->[1], $pos->[2], $pos->[3]->ipot, $tag, sqrt($pos->[4]) ); $seen{$pos->[3]->tag}->[1] = sqrt($pos->[4]); }; return $string; }; sub xyz_list { my ($self, $pattern) = @_; $pattern ||= " %-10s %9.5f %9.5f %9.5f\n"; my $string = q{}; my @list = @ {$self->cluster }; foreach my $pos (@list) { $string .= sprintf($pattern, $pos->[3]->element, $pos->[0], $pos->[1], $pos->[2]); }; return $string; }; sub alchemy_list { my ($self, $pattern) = @_; $pattern ||= " %4d %-2s %9.5f %9.5f %9.5f 0.00\n"; my $string = q{}; my @list = @ {$self->cluster }; my $count = 1; foreach my $pos (@list) { $string .= sprintf($pattern, $count, $pos->[3]->element, $pos->[0], $pos->[1], $pos->[2]); ++$count; }; return $string; }; sub R { my ($self, $x, $y, $z) = @_; return sqrt($x**2 + $y**2 + $z**2); }; sub potentials_list { my ($self, $pattern) = @_; $self->set_ipots if (not $self->is_ipots_set); $pattern ||= " %d %-2d %-10s\n"; my $is_feff8 = 0; my @list = split(" ", $pattern); $is_feff8 = 2 if ($#list > 3); my ($cell, $core) = $self->get("cell", "core"); my @sites = @{ $cell->sites }; my $string = q{}; my %seen = (); my ($abs) = $cell->central($core); my $l = Xray::Absorption->get_l($abs->element); @list = ($is_feff8) ? (0, get_Z($abs->element), $abs->element, $l, $l, 0.001) : (0, get_Z($abs->element), $abs->element); $string .= sprintf($pattern, @list); foreach my $s (sort {$a->ipot <=> $b->ipot} @sites) { next if not $s->ipot; next if $seen{$s->ipot}; $l = Xray::Absorption->get_l($s->element); @list = ($is_feff8) ? ($s->ipot, get_Z($s->element), $s->element, $l, $l, $s->stoi) : ($s->ipot, get_Z($s->element), $s->element); $string .= sprintf($pattern, @list); $seen{$s->ipot} = 1; }; return $string; }; sub sites_list { my ($self, $rhash) = @_; $self->populate if (not $self->is_populated); my $prec = '%'.$self->co->default("atoms", "precision"); $rhash->{pattern} ||= " %-2s $prec $prec $prec %-10s\n"; $rhash->{prefix} ||= q{}; my $cell = $self->cell; my $rlist = $cell->sites; my $string = q{}; foreach my $l (@$rlist) { $string .= $rhash->{prefix} . sprintf($rhash->{pattern}, ucfirst(lc($l->element)), $l->x, $l->y, $l->z, $l->tag); }; return $string; }; sub p1_list { my ($self, $rhash) = @_; $self->populate if (not $self->is_populated); my $prec = '%'.$self->co->default("atoms", "precision"); $rhash->{pattern} ||= " %-2s $prec $prec $prec %-10s\n"; $rhash->{prefix} ||= q{}; $rhash->{cartesian} ||= 0; my $cell = $self->cell; my $rlist = $cell->contents; my $string = q{}; foreach my $l (@$rlist) { if ($rhash->{cartesian}) { my @list = $cell->metric($$l[1], $$l[2], $$l[3]); $string .= $rhash->{prefix} . sprintf($rhash->{pattern}, ucfirst(lc($l->[0]->element)), @list, $l->[0]->tag) } else { $string .= $rhash->{prefix} . sprintf($rhash->{pattern}, ucfirst(lc($l->[0]->element)), $$l[1], $$l[2], $$l[3], $l->[0]->tag) }; }; return $string; }; sub sg { my ($self, $which, $pattern) = @_; $self->populate if (not $self->is_populated); ($which = "shorthand") if ($which eq 'nicknames'); my $cell = $self->cell; my $rhash = $cell->group->data; $pattern ||= " %-8s %-8s %-8s\n"; my ($prefix, $postfix) = ($which =~ m{(?:bravais|shiftvec)}) ? (" ", $/) : (q{}, q{}); ## typo? return q{} if (not is_SpaceGroup($which)); ## --------------------------------- ## number of positions if ($which eq "npos") { my @positions; if (defined($$rhash{positions})) { @positions = @ {$$rhash{positions}}; } else { @positions = @ {$$rhash{$cell->group->{setting}}}; } return $#positions + 1; }; ## --------------------------------- ## number or symbol #return $$rhash{$which} if ($which =~ m{(?:number|full|newsymbol|thirtyfive)}); if (($which =~ m{(?:number|newsymbol|thirtyfive)})) { return ucfirst($self->cell->group->$which) || ''; }; if ($which eq 'full') { return ucfirst($self->cell->group->full); }; ## --------------------------------- ## key is absent from this entry in database return "$prefix$postfix" if ((not exists($$rhash{$which})) and ($which ne "bravais") and ($which ne "positions")); ## --------------------------------- ## schoenflies return ucfirst($$rhash{schoenflies}) if ($which eq "schoenflies"); ## --------------------------------- ## nicknames return join(", ", @{$$rhash{shorthand}}) if ($which eq "shorthand"); ## --------------------------------- ## shift vector from Int'l Tables if ($which eq "shiftvec") { #my @shift = map {fract($FRAC*$_, $FRAC)} @{ $$rhash{shiftvec} }; #return sprintf($pattern, map {$_->as_mixed_string} @shift); my @shift = map {$self->fract($_)} @{ $$rhash{shiftvec} }; return sprintf($pattern, @shift); }; ## --------------------------------- ## Bravais translations if ($which eq "bravais") { my @bravais = @{ $cell->group->bravais }; my $string = q{}; while (@bravais) { #my @vec = (fract($FRAC*shift(@bravais), $FRAC), # fract($FRAC*shift(@bravais), $FRAC), # fract($FRAC*shift(@bravais), $FRAC), # ); #$string .= sprintf($pattern, map {$_->as_mixed_string} @vec); my @vec = ($self->fract(shift(@bravais)), $self->fract(shift(@bravais)), $self->fract(shift(@bravais)), ); $string .= sprintf($pattern, @vec); }; return $string; }; ## --------------------------------- ## symetric positions if ($which eq "positions") { my @positions; if (defined($$rhash{positions})) { @positions = @ {$$rhash{positions}}; } else { @positions = @ {$$rhash{$cell->group->{setting}}}; } my $string = q{}; my $npos = $#positions + 1; #$string .= " $npos positions:\n"; foreach my $pos (@positions) { my @this = @{ $pos }; map { $this[$_] =~ s{\$}{}g } (0 .. 2); $string .= sprintf($pattern, map {($_ =~ m{\A\-}) ? $_ : " $_"} @this); }; return $string; }; return q{}; }; sub all_titles { my ($self, $prefix) = @_; $prefix ||= " TITLE "; my @titles = @{ $self->titles }; my $string = q{}; foreach my $t (@titles) { $string .= $prefix . $t . $/; }; return $string; }; sub update_absorption { my ($self) = @_; $self->_absorption if not $self->absorption_done; $self->_mcmaster if not $self->mcmaster_done; $self->_i0 if (($self->gases_set) and not $self->i0_done); $self->_self if (($self->gases_set) and not $self->self_done); return $self; }; sub update_edge { my ($self) = @_; my $edge = $self->edge; my ($central, $xcenter, $ycenter, $zcenter) = $self -> cell -> central($self->core); ##print $self->core, $/; ##print $central, $/; ##print join(" ", $central->meta->get_attribute_list), $/; #Demeter->trace; if ($edge) { $self->edge($edge); # be sure to trigger assignment of eedge } else { my $z = get_Z( $central->element ); ($z > 57) ? $self->edge('l3') : $self->edge('k'); }; return $self; }; sub handle_colon { my ($self) = @_; my $g = $self->space; if ($self->is_first) { ## remove :1 from the space group symbol $g =~ s{:[12]}{}; $self->space($g); ## set the shift vector according to the entry in the space group database $self->shiftvec($self->cell->group->shiftvec); } elsif ($self->is_second) { ## remove :2 from the space group symbol $g =~ s{:[12]}{}; $self->space($g); }; return $self; }; sub Write { my ($self, $type) = @_; $type ||= "feff6"; $type = lc($type); ($type = 'feff6') if ($type eq 'feff'); $self->update_absorption; return $self->atoms_file if ($type eq 'atoms'); return $self->atoms_file('p1') if ($type eq 'p1'); return $self->template('overfull') if ($type eq 'overfull'); return $self->template("absorption") if (($type eq 'absorption') and $self->gases_set); return $self->template("mcmaster") if (($type eq 'absorption') and not $self->gases_set); if ($type eq 'spacegroup') { $self->populate if (not $self->is_populated); return $self->spacegroup_file(0, '# '); }; return $self->Write_feff($type) if ($type =~ m{feff|xyz|alc}); ## still need: overfull, p1_cartesian, gnxas ## fallback return $self->Write_feff('feff6'); }; sub Write_feff { my ($self, $type) = @_; $self->build_cluster if (not $self->is_expanded); return q{} if ($#{$self->cluster} == -1); return q{} if ($#{$self->cluster} == 0); my $string = q{}; if ($type =~ m{feff}) { # none of this header stuff for xyz or alchemy or the like $string .= $self->template('copyright', {type=> $type, prefix => ' * '}); if ($self->co->default("atoms", "atoms_in_feff")) { $string .= $self->template('prettyline', {prefix => ' * '}); $string .= $self->atoms_file('feff', ' * '); $string .= $self->template('prettyline', {prefix => ' * '}); $string .= $/; } if ($self->gases_set) { $string .= $self->template('absorption', {prefix => ' * '}); } else { $string .= $self->template('mcmaster', {prefix => ' * '}); }; }; $string .= $self->template($type); return $string; }; sub atoms_file { my ($self, $is_p1, $prefix) = @_; $is_p1 ||= 0; $prefix ||= q{}; $self->populate if (not $self->is_populated); my $cell = $self -> cell; my $string = q{}; my $type = ($is_p1 eq 'p1') ? 'P1' : ($is_p1) ? q{} : 'Atoms'; $string .= $self->template("copyright", {prefix=>$prefix, type=>$type}) if $type; $string .= $self->template("atoms_header", {prefix=>$prefix, is_p1=>($type eq 'P1')}); $string .= ($type eq 'P1') ? $self->p1_list({prefix=>$prefix}) : $self->sites_list({prefix=>$prefix}); return $string; }; sub overfull_list { my ($self, $nclus) = @_; $self->populate if (not $self->is_populated); my $cell = $self -> cell; my $p1cell = $self->p1_list(); my $string = q{}; my $count = 0; foreach my $line (split($/, $p1cell)) { my ($elem, $x, $y, $z, $tag) = split(" ", $line); $string .= sprintf(" %-2s %9.5f %9.5f %9.5f\n", $elem, $cell->metric($x, $y, $z)); ++$count; ## atoms with a near-0 coordinate if (abs($x) < Demeter->co->default('atoms', 'overfull_margin')) { $string .= sprintf(" %-2s %9.5f %9.5f %9.5f\n", $elem, $cell->metric($x+1, $y, $z)); ++$count; }; if (abs($y) < Demeter->co->default('atoms', 'overfull_margin')) { $string .= sprintf(" %-2s %9.5f %9.5f %9.5f\n", $elem, $cell->metric($x, $y+1, $z)); ++$count; }; if (abs($z) < Demeter->co->default('atoms', 'overfull_margin')) { $string .= sprintf(" %-2s %9.5f %9.5f %9.5f\n", $elem, $cell->metric($x, $y, $z+1)); ++$count; }; if ((abs($x) < Demeter->co->default('atoms', 'overfull_margin')) and (abs($y) < Demeter->co->default('atoms', 'overfull_margin'))) { $string .= sprintf(" %-2s %9.5f %9.5f %9.5f\n", $elem, $cell->metric($x+1, $y+1, $z)); ++$count; }; if ((abs($x) < Demeter->co->default('atoms', 'overfull_margin')) and (abs($z) < Demeter->co->default('atoms', 'overfull_margin'))) { $string .= sprintf(" %-2s %9.5f %9.5f %9.5f\n", $elem, $cell->metric($x+1, $y, $z+1)); ++$count; }; if ((abs($y) < Demeter->co->default('atoms', 'overfull_margin')) and (abs($z) < Demeter->co->default('atoms', 'overfull_margin'))) { $string .= sprintf(" %-2s %9.5f %9.5f %9.5f\n", $elem, $cell->metric($x, $y+1, $z+1)); ++$count; }; if ((abs($x) < Demeter->co->default('atoms', 'overfull_margin')) and (abs($y) < Demeter->co->default('atoms', 'overfull_margin')) and (abs($z) < Demeter->co->default('atoms', 'overfull_margin'))) { $string .= sprintf(" %-2s %9.5f %9.5f %9.5f\n", $elem, $cell->metric($x+1, $y+1, $z+1)); ++$count; }; ## atoms with a near-1 coordinate if (abs(1-$x) < Demeter->co->default('atoms', 'overfull_margin')) { $string .= sprintf(" %-2s %9.5f %9.5f %9.5f\n", $elem, $cell->metric($x-1, $y, $z)); ++$count; }; if (abs(1-$y) < Demeter->co->default('atoms', 'overfull_margin')) { $string .= sprintf(" %-2s %9.5f %9.5f %9.5f\n", $elem, $cell->metric($x, $y-1, $z)); ++$count; }; if (abs(1-$z) < Demeter->co->default('atoms', 'overfull_margin')) { $string .= sprintf(" %-2s %9.5f %9.5f %9.5f\n", $elem, $cell->metric($x, $y, $z-1)); ++$count; }; if ((abs(1-$x) < Demeter->co->default('atoms', 'overfull_margin')) and (abs(1-$y) < Demeter->co->default('atoms', 'overfull_margin'))) { $string .= sprintf(" %-2s %9.5f %9.5f %9.5f\n", $elem, $cell->metric($x-1, $y-1, $z)); ++$count; }; if ((abs(1-$x) < Demeter->co->default('atoms', 'overfull_margin')) and (abs(1-$z) < Demeter->co->default('atoms', 'overfull_margin'))) { $string .= sprintf(" %-2s %9.5f %9.5f %9.5f\n", $elem, $cell->metric($x-1, $y, $z-1)); ++$count; }; if ((abs(1-$y) < Demeter->co->default('atoms', 'overfull_margin')) and (abs(1-$z) < Demeter->co->default('atoms', 'overfull_margin'))) { $string .= sprintf(" %-2s %9.5f %9.5f %9.5f\n", $elem, $cell->metric($x, $y-1, $z-1)); ++$count; }; if ((abs(1-$x) < Demeter->co->default('atoms', 'overfull_margin')) and (abs(1-$y) < Demeter->co->default('atoms', 'overfull_margin')) and (abs(1-$z) < Demeter->co->default('atoms', 'overfull_margin'))) { $string .= sprintf(" %-2s %9.5f %9.5f %9.5f\n", $elem, $cell->metric($x-1, $y-1, $z-1)); ++$count; }; }; $self->nclus($count); $string = $count . $/ . $string if ($nclus); return $string; }; sub spacegroup_file { my ($self) = @_; my $prefix = '# '; my $string = $self->template("copyright", {prefix=>$prefix, type=>"space group"}); $string .= $self->template("atoms_header", {prefix=>$prefix}); $string .= $self->sites_list({prefix=>$prefix}); $string .= $/; $string .= $self->template("spacegroup"); return $string; }; sub is_polarization { my ($self, $token) = @_; return '' if (($self->polarization->[0]>$EPSILON4) or ($self->polarization->[1]>$EPSILON4) or ($self->polarization->[2]>$EPSILON4)); return $token || ''; }; override serialization => sub { my ($self) = @_; my %cards = (); foreach my $key (qw(space a b c alpha beta gamma rmax rpath rss edge iedge eedge core corel partial_occupancy shiftvec cif record titles ipot_style nitrogen argon krypton xenon helium gases_set xsec deltamu density mcmaster i0 selfamp selfsig netsig is_imported is_populated is_ipots_set is_expanded is_rhomb is_hex is_first is_second absorption_done mcmaster_done i0_done self_done nclus polarization ellipticity sites)) { # sites cluster $cards{$key} = $self->$key; }; my $text = YAML::Tiny::Dump(\%cards); return $text; }; __PACKAGE__->meta->make_immutable; 1; =head1 NAME Demeter::Atoms - Convert crystallographic data to atomic lists =head1 VERSION This documentation refers to Demeter version 0.9.26. =head1 SYNOPSIS These lines behave much like any other version of Atoms: use Demeter; my $atoms = Demeter::Atoms -> new() -> read_inp($ARGV[0]||"atoms.inp") -> Write('feff6'); =head1 DESCRIPTION This module implements Atoms in the Demeter system. The purpose of Atoms is to convert crystallographic data into a list of atomic coordinates of the sort used by Feff as input data. This greatly simplifies the chore of making Feff input files for crystalline materials. =head1 ATTRIBUTES The following are the attributes of the Data object. Attempting to access an attribute not on this list will throw an exception. The type of argument expected in given in parentheses. i.e. number, integer, string, and so on. The default value, if one exists, is given in square brackets. =head2 Input parameters =over 4 =item C (string) The space group ofthe crystal. This can be in any form recognized by L, including Hermann-Maguin, Schoenflies, number, or one of a few nicknames. =item C (real) The length of the A lattice constant in Angstroms. =item C (real) The length of the B lattice constant in Angstroms. =item C (real) The length of the C lattice constant in Angstroms. =item C (real) The angle between B and C, in degrees. =item C (real) The angle between A and C, in degrees. =item C (real) The angle between A and B, in degrees. =item C (real) The extent of the cluster generated from the input crystal data. =item C (real) The value used for the RMAX keyword in the F file. This is the length of the longest path to be calculated by the pathfinder. A value much larger than about 6 will bog down Demeter's pathfinder in its current form. =item C (string) [k or l3 depending on Z-number] The edge of the absorber. =item C (string) [first sitin the list] The identifier of the absorber. This should be one of the site tags. =item C (vector) [0,0,0] The value of the shift vector, should one be necessary. =item C (filename) The name of an atoms input file. =item C (filename) The name of a CIF file. =item C (string) [0] The record to import from a multi-record CIF file. The default is to read the first record. Note that this is zero-based while you user interface probably should be one-based. =item C (array of strings) An array of strings containing the title lines read from the input data. =item C