2
######################################################################
3
## This is the Xray::Atoms.pm module. It contains various subroutines
4
## used by the various versions of atoms. See the pod for details and
5
## one of the atoms versions for how they are implemented.
7
## This program is copyright (c) 1998-2006 Bruce Ravel
9
## http://cars9.uchicago.edu/~ravel/software/
11
## -------------------------------------------------------------------
12
## All rights reserved. This program is free software; you can
13
## redistribute it and/or modify it under the same terms as Perl
16
## This program is distributed in the hope that it will be useful,
17
## but WITHOUT ANY WARRANTY; without even the implied warranty of
18
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19
## Artistic License for more details.
20
## -------------------------------------------------------------------
22
######################################################################
23
## Time-stamp: <06/01/02 11:56:22 bruce>
24
######################################################################
29
Xray::Atoms - Utilities and data structures for the Atoms program
34
use Xray::Atoms qw(parse_input keyword_defaults parse_atp
39
This module contains the utility subroutines used by the program
40
Atoms. It defines a number of exportable routines for performing
41
chores specific to Atoms. It also defines an object for storing input
42
data to the various programs which use this file. For details about
43
the crystallographic objects used by Atoms, see L<Xray::Xtal>.
45
There are several versions of Atoms with different user interfaces all
46
of which use routines from this module. These include the command
47
line, Tk, and CGI versions.
55
use vars qw($VERSION $atoms_dir $languages @available_languages
56
$cvs_info $module_version $messages @ISA @EXPORT @EXPORT_OK);
60
@ISA = qw(Exporter AutoLoader);
61
# Items to export into callers namespace by default. Note: do not export
62
# names by default without a very good reason. Use EXPORT_OK instead.
63
# Do not simply export all your public functions/methods/constants.
64
@EXPORT_OK = qw(build_cluster number rcdirectory rcfile_name
65
absorption mcmaster i_zero self);
67
$VERSION = '3.0beta10';
68
$cvs_info = '$Id: Atoms.pm,v 1.40 2001/12/21 02:28:44 bruce Exp $ ';
69
$module_version = 1.41; #(split(' ', $cvs_info))[2];
75
use Chemistry::Elements qw(get_name get_Z get_symbol);
76
use Statistics::Descriptive;
79
use Ifeffit::FindFile;
80
my $STAR_Parser_exists = (eval "require STAR::Parser");
81
if ($STAR_Parser_exists) {
83
require STAR::DataBlock;
84
import STAR::DataBlock;
87
#use constant EV2RYD => 13.605698;
88
use constant EPSILON => 0.00001;
90
my $elem_match = '([BCFHIKNOPSUVWY]|A[cglrstu]|B[aeir]|C[adelorsu]|Dy|E[ru]|F[er]|G[ade]|H[efgo]|I[nr]|Kr|L[aiu]|M[gno]|N[abdeip]|Os|P[abdmortu]|R[abehnu]|S[bceimnr]|T[abcehilm]|Xe|Yb|Z[nr])';
92
## location of atp files
94
use vars qw($atp_dir $lib_dir $is_windows);
95
$is_windows = $Ifeffit::FindFile::is_windows;
96
$atoms_dir = Xray::Xtal::identify_self();
97
$atp_dir = ($is_windows) ?
98
Ifeffit::FindFile->find("atoms", "atp_sys") :
99
File::Spec->catfile($atoms_dir, "atp");
100
$lib_dir = ($is_windows) ?
101
Ifeffit::FindFile->find("atoms", "xray_lib") :
102
File::Spec->catfile($atoms_dir, "lib");
104
## my $languagerc = (($^O eq 'MSWin32') or ($^O eq 'cygwin')) ?
105
## File::Spec->catfile($lib_dir, "languages");
106
## eval "do '$languagerc'" or warn "Language rc file not found. Using English.$/";
109
use vars qw(%meta %colors %fonts);
110
## use vars qw($always_write_feff $atoms_language $write_to_pwd
111
## $prefer_feff_eight $absorption_tables $dafs_default
112
## $plotting_hook $default_filepath $unused_modifier
113
## $display_balloons $no_crystal_warnings $one_frame $convolve_dafs
114
## $never_ask_to_save $ADB_location);
115
%meta = (ADB_location => "http://cars9.uchicago.edu/atomsdb/",
116
absorption_tables => 'elam',
117
always_write_feff => 0,
118
atoms_language => 'english',
120
dafs_default => 'cl',
121
default_filepath => '',
122
display_balloons => 1,
123
never_ask_to_save => 0,
124
no_crystal_warnings => 0,
127
prefer_feff_eight => 0,
128
unused_modifier => 'Shift',
132
## keep these refs to anon hashes so as not to break rcfiles from
133
## alpha16 and earlier
134
## use vars qw(); # unused variables from atomsrc
135
## ## must declare variables for colors and fonts for TkAtoms
136
## use vars qw($c_foreground $c_background $c_trough $c_entry $c_label
137
## $c_balloon $c_button $c_buttonActive $c_sgbActive $c_sgbGroup
138
## $c_done $c_todo $c_plot
139
## $f_balloon $f_label $f_menu $f_button $f_header $f_entry $f_sgb);
141
use vars qw($rcfile);
142
$rcfile = (($^O eq 'MSWin32') or ($^O eq 'cygwin')) ?
143
File::Spec->catfile($lib_dir, "tkatoms.ini"):
144
File::Spec->catfile($lib_dir, "atomsrc");
145
(-e $rcfile) and &read_rc($rcfile);
146
my $users_rc = &rcfile_name;
147
(-e $users_rc) and &read_rc($users_rc);
148
$meta{atoms_language} = lc($meta{atoms_language});
149
$meta{absorption_tables} = lc($meta{absorption_tables});
150
Xray::Absorption -> load($meta{absorption_tables});
151
#Xray::Absorption -> load('Elam');
152
$meta{unused_modifier} = ucfirst($meta{unused_modifier});
153
($meta{unused_modifier} =~ /Alt|Control|Meta|Shift/) or
154
$meta{unused_modifier} = 'Shift';
157
if (not $languages) {
158
$language_file = (($^O eq 'MSWin32') or ($^O eq 'cygwin')) ?
159
File::Spec->catfile($lib_dir, "atomsrc.en"):
160
File::Spec->catfile($lib_dir, "atomsrc.en");
162
$language_file = "atomsrc." . $$languages{$meta{atoms_language}};
163
$language_file = (($^O eq 'MSWin32') or ($^O eq 'cygwin')) ?
164
File::Spec->catfile($lib_dir, $language_file):
165
File::Spec->catfile($lib_dir, $language_file);
166
unless (-e $language_file) {
167
$language_file = (($^O eq 'MSWin32') or ($^O eq 'cygwin')) ?
168
File::Spec->catfile($lib_dir, "atomsrc.en"):
169
File::Spec->catfile($lib_dir, "atomsrc.en");
172
eval "do '$language_file'"; # read strings
173
($meta{atoms_language} eq 'english') or
174
Xray::Xtal::set_language($meta{atoms_language});
176
my $overfull_margin = 0.1;
178
# Preloaded methods go here.
180
=head1 THE KEYWORDS OBJECT
182
To simplify the handling of input data, this module provides an
183
object. It is used like this:
185
my $keywords = Xray::Atoms -> new();
186
$keywords -> make('rmax' => 5.7)
188
The make method is the general way of storing data in the keywords
189
object. In general, this method assumes that the keyword takes a
190
scalar value. The data structure is just a hash, so the keyword can
191
be anything. Indeed, keyword recognition is not done in this module.
192
Keyword/value pairs are just blindly stored for later processing.
193
Several keywords are treated specially because they are commonly used
194
in the programs included in the Atoms package. These keywords
196
a b c alpha beta gamma argon krypton nitrogen rmax
198
are treated as numbers. This means that they are eval-ed (and so can
199
be short expressions like 1/2 or 3.5+0.001) and return 0 if they
200
cannot be interpreted as numbers. This evaluation is done in a
203
A couple of keywords are treated specially.
205
$keywords -> make(title=>'This is a title!');
207
pushes the given string onto an anonymous array of strings. In this
208
way, any number of title lines can be kept in the object.
210
The C<shift> keyword is also handled specially. This takes three
211
arguments, the three components of the shift vector. A syntax like
214
$keywords -> make('shift' => $x, $y, $z);
216
This stores the three values of the shift vector as an anonymous
224
my $classname = shift;
227
$self->{'title'} = [];
228
$self->{'edge'} = "";
229
$self->{'core'} = "";
230
$self->{'argon'} = 0;
231
$self->{'krypton'} = 0;
232
$self->{'nitrogen'} = 0;
234
$self->{'shift'} = [0,0,0];
236
$self->{'space'} = '';
240
$self->{'alpha'} = 0;
242
$self->{'gamma'} = 0;
244
$self->{'sites'} = [];
245
# operational keywords
246
$self->{'overfull_margin'} = 0.1;
247
$self->{'found_output'} = 0;
248
$self->{'quiet'} = 0;
249
$self->{'program'} = 'atoms';
252
my @atp_dir_list = ();
253
push @atp_dir_list, $Xray::Atoms::atp_dir;
254
push @atp_dir_list, File::Spec->catfile(&rcdirectory, "atp");
255
foreach my $dir (@atp_dir_list) {
257
opendir (ATPDIR, $dir) ||
258
die $$messages{'no_directory'} . " " . $dir . $/;
259
push @atpfiles, grep /\.atp/, readdir ATPDIR;
262
}; ## feff feff8 p1 unit alchemy xyz geom symmetry test
263
$self->{'files'} = {};
265
$self->{'ipots'} = 'species'; # 'sites'
266
$self->{'ipot_vals'} = {};
267
$self->{'always_feff'} = $meta{always_write_feff};
268
$self->{'prefer_feff_8'} = $meta{prefer_feff_eight};
269
$self->{'language'} = $meta{atoms_language};
270
$self->{'write_to_pwd'} = $meta{write_to_pwd};
272
bless($self, $classname);
276
##"fdat" => 0, "refile" => "reflect.dat",
277
##"qvec" => (0,0,0), "nepoints" => 100,
278
##"feout" => "fe.dat", "egrid" => 1,
279
##"reflections" => (0,0,0), "noanomalous" => 0 );
280
## corrections, index
285
unless (($#_ % 2) || (grep /\b(files|sites)\b/, @_)) {
286
my $this = (caller(0))[3];
287
croak "$this " . $$messages{make_error};
291
my $die = $self->{die} || 1;
296
($att eq 'title') && do {
297
foreach my $l (split(/\n/, $value)) {
298
push @{$self->{'title'}}, $l unless ($l =~ /^\s*$/);
302
($att =~ /(shift|qvec)/) && do {
303
my $v1 = number($value, $die, $self);
304
my $v2 = number(shift, $die, $self);
305
my $v3 = number(shift, $die, $self);
306
$self->{$1} = [$v1, $v2, $v3];
309
($att =~ /^(argon|krypton|nitrogen|rmax)$/) && do {
310
$self->{$1} = number($value, $die, $self);
313
($att =~ /^(a|b|c|alpha|beta|gamma)$/) && do {
314
$self->{$1} = number($value, $die, $self);
317
($att eq 'atp') && do {
319
$self->{atp}{$value} = 1;
322
($att eq 'files') && do {
325
($v2 =~ /^\s*$/) and $v2 = undef;
326
$self->{files}{$value} = $v2;
329
($att eq 'sites') && do {
331
my $x = number(shift, $die, $self);
332
my $y = number(shift, $die, $self);
333
my $z = number(shift, $die, $self);
335
my $o = number(shift, $die, $self);
336
push @{$self->{'sites'}}, [$e, $x, $y, $z, $t, $o];
340
$self->{$att} = $value;
354
=head2 C<parse_input>
356
This method is used by atoms to read an atoms input file. It is
359
$keywords -> parse_input($file, 0);
361
The first argument is the name of the input file. The input file is
362
presumed to be of the atoms.inp sort unless the file extension is
363
C<.cif>, in which case it is presumed to be a CIF file. This check is
364
made case-insensitively.
366
The second argument is tells the error trapping mechnism what kind of
367
environment you are working in. If you have written a command-line
368
program, it should be 0. In a Tk program, it should 1. In a CGI
369
script it should be 2.
371
The contents of the file are parsed and stored in a keywords object.
372
The method is fairly clever about interpreting the file and offering
373
warning and error messages. It interprets numbers in a fairly secure
374
manner, thus allowing simple math expressions to be used in the atom
375
list and as the shift vectors.
377
Alternately, you can pass a CIF file with the syntax
379
$keywords -> parse_input($file, 0, 'cif');
383
$keywords -> parse_input($file, 0, 'cif', 2);
385
The third argument is "cif" or "inp" to specify the input file type.
386
The fourth argument tells the CIF parser which structure to use from a
387
multi-structure file.
391
$$::messages{adb_unknown} = "You requested an unknown ADB file";
394
my $href = abbrev qw(a alpha argon atoms b basis beta c core corrections
395
edge egrid emax emin egrid estep dopants fdat feff
396
feff8 gamma geom index ipots krypton
397
nepoints nitrogen noanomalous out output p1
398
qvec refile reflection reflections rmax shift
399
space title unit xanes);
402
my ($file, $die, $type, $entry) = @_;
404
$type = ($type =~ /cif/i) ? "cif" : "inp";
406
my ($nsites, $ntitles, $iline) = (0,0,0);
408
## is this a CIF file?
409
my ($nn,$pp,$suff) = fileparse($file,
410
".inp", ".INP", ".Inp",
411
".cif", ".CIF", ".Cif");
412
cif($keys, $file ,$entry), return
413
if ($STAR_Parser_exists and ((lc($type) eq 'cif') or ($suff =~ /cif$/i)));
416
## divy up the keywords into birds of
418
my @normal = ("argon", "center", "core", "ipots", "krypton",
421
my @threevecs = ("shift", "dafs", "qvec", "reflection");
422
my @cellwords = ("a", "alpha", "b", "beta", "c", "gamma");
423
my @output_types = ("feff", "feff8", "p1", "geom", "unit", "out");
424
my @dafs = ("emin", "emax", "estep", "egrid");
425
my @deprecated = ("fdat", "nepoints", "xanes", "modules",
426
"message", "noanomalous", "self", "i0", "mcmaster",
427
"dwarf", "reflections", "refile",
428
"egrid", "index", "corrections");
429
my ($term, $prompt, $stdin);
431
if ($file =~ /^http:/) {
432
if (eval "require LWP::Simple") {
433
unless (LWP::Simple::head($file)) { # handle unknown ADB file
434
$keys -> warn_or_die("$$::messages{adb_unknown}:\n\t$file.\n", $die);
437
open ($fh, "GET $file |")
438
or die "A pipe using LWP::Simple could not be opened to fetch ADB file.\n";
440
die "You must install LWP::Simple to fetch ADB files.\n";
442
} elsif ($file ne '____stdin') {
443
open ($fh, $file) || die "could not open " . $file . " for reading" . $/;
445
## what about STDIN???
447
READ: while (<$fh>) {
451
next if /^\s*$/; # skip blanks
452
next if /^\s*[!\#%*]/; # skip comment lines
454
(m/<HTML>/) and do { # handle an unknown ADB file name
457
(my $line = $_) =~ s/^\s+//; # trim leading blanks
458
$line =~ s/\r//g; # strip spurious control-M characters
459
@line = split(/\s*[ \t=,]\s*/, $line);
461
KEYWORDS: while (@line) {
462
my $word = lc(shift @line);
463
next READ if ($word =~ /^[!\#%*]/); # skip comment lines
465
$word = $href -> {$word};
468
## my $message = "\"$wasword\": " . $$messages{'unknown_keyword'} .
469
## ", " . $$messages{'line'} . " $..$/";
470
## $keys -> warn_or_die($message, $die);
473
($word eq "end") && last READ;
477
($word =~ /^(com|tit)/) && do {
478
## trim `title' or `comment' and leading spaces
480
my $string = substr($line, length($word));
481
$string =~ s/^\s+//; $string =~ s/^=//; $string =~ s/^\s+//;
482
unless ($keys -> {"quiet"}) {
483
print STDOUT " ", $$messages{'title'}, " > ",
484
substr($string, 0, 59), "$/";
486
$keys -> make('title'=>$string);
487
foreach my $a (@line) {
492
## edge/hole can be error-checked en
494
($word =~ /^(edg|hol)/) && do {
495
my $value = shift @line;
496
$keys -> make('edge'=>$value);
499
## keywords which take the next word
501
## argon center core krypton nitrogen rmax
502
(grep(/\b$word\b/, @normal)) && do {
503
($word eq "center") && ($word = "core");
504
my $value = shift @line;
505
$keys -> make($word=>$value);
508
## DAFS keywords which take the next
509
## word literally: emin estep emax
510
(grep(/\b$word\b/, @dafs)) && do {
511
($word eq "egrid") && ($word = "estep");
512
my $value = shift @line;
513
$keys -> make($word=>$value);
516
## three vectors, take the next three
518
(grep(/\b$word\b/, @threevecs)) && do {
520
$value[0] = shift @line;
521
$value[1] = shift @line;
522
$value[2] = shift @line;
523
($word eq "dafs") && ($word = "qvec");
524
($word eq "reflection") && ($word = "qvec");
526
## ----- evaluate the coordinates safely
527
## see `perldoc Safe' for details
528
## my $message = $$messages{'tainted_vector'} . $/;
529
## $message .= $$messages{'file'} . ": $file, ";
530
## $message .= $$messages{'line'} . ": $.$/";
531
## my $cpt = new Safe;
532
## ##$cpt->share('@value');
533
## $ {$cpt->varglob('x')} = $cpt->reval($value[0]);
534
## $ {$cpt->varglob('y')} = $cpt->reval($value[1]);
535
## $ {$cpt->varglob('z')} = $cpt->reval($value[2]);
536
## @value = ($ {$cpt->varglob('x')}, $ {$cpt->varglob('y')},
537
## $ {$cpt->varglob('z')} );
538
## test_safe_return($message, @value);
540
## this is a much weaker form of security, but it works with
542
@value = map {number($_,$die,$keys)} @value;
543
## ----- end of safe evaluation
545
$keys -> make($word=>$value[0],$value[1],$value[2]);
548
## logical flags and diagnostics take
549
## boolean values, on if next word is
551
##(grep(/\b$word\b/, @logicals)) && do {
552
## my $value = shift @line;
553
## $keywords{$word} = ($value =~ /^(t|y|on)/) ? 1 : 0;
556
## output styles take filenames
557
($word =~ /^outp/) && do {
558
my $type = shift @line;
559
my $value = shift @line;
560
#push @$outputs, $type;
561
if ($keys->{'program'} eq 'dafs') {
562
undef $keys->{atp}; undef $keys->{files};
564
$keys -> make('atp' => $type);
565
$keys -> make('files' => $type, $value);
566
$keys -> make("found_output" => 1);
569
## backwards compatability
570
(grep(/\b$word\b/, @output_types)) && do {
571
my $value = shift @line;
572
if ($word eq "out") {
573
$keys -> make('atp'=>'feff', 'feff'=>$value);
574
#push @$outputs, "feff";
576
my $set = ($value =~ /^(t|y|on)/);
578
$keys -> make('atp'=>$word);
579
## ($word eq "feff") && ($keys -> make('feff' => "feff.inp"));
580
## ($word eq "feff8") && ($keys -> make('feff8'=> "feff.inp"));
581
## ($word eq "geom") && ($keys -> make('geom' => "geom.dat"));
582
## ($word eq "p1") && ($keys -> make('p1' => "p1.inp" ));
583
## ($word eq "unit") && ($keys -> make('unit' => "unit.dat"));
584
#push @$outputs, $word;
589
## deprecated keywords
590
(grep(/\b$word\b/, @deprecated)) && do {
592
my $mess = "\"$word\": " . $$messages{'deprecated'} . ", " .
593
$$messages{'line'} . " $..$/";
595
($word =~ /dafs|reflections|qvec/) && do {
596
shift @line; shift @line;
600
## space group symbol
601
($word =~ /^spa/) && do {
602
my $lline = lc($line);
603
my $space = substr($line, index($lline,"space")+6);
604
$space =~ s/^[\s=,]+//;
605
$space = substr($space, 0, 10); # next 10 characters
606
$space =~ s/[!\#%*].*$//; # trim off comments
607
$keys -> make('space'=>$space);
608
foreach my $x (split " ", $space) {
613
## a,b,c,alpha,beta,gamma
614
(grep(/\b$word\b/, @cellwords)) && do {
615
my $value = shift @line;
616
$value = sprintf "%11.7f", $value;
617
my $attribute = lc $word;
618
$keys -> make($attribute=>$value);
619
## the angle attribute carries the most recently set angle
620
## this is needed to correctly interpret monoclinic cells
621
## with all angles equal to 90 (pathology!!)
622
##($attribute =~ /(alpha|beta|gamma)/i) &&
623
## ($cell -> make( "Angle"=>$1 ) );
627
($word =~ /^dop/) && do {
628
my $mess = $$messages{'dopants'} . ", " . $$messages{'line'} . " $..$/";
629
$keys -> warn_or_die($mess, $die);
630
shift @line; shift @line; shift @line;
633
($word =~ /^bas/) && do {
634
croak $$messages{'basis'}, $/;
637
($word =~ /^ato/) && do {
641
(! $_) && last KEYWORDS;
643
last KEYWORDS if /---/;
644
next ATOMS if /^\s*$/; # skip blanks
645
next ATOMS if /^\s*[!\#%*]/; # skip comment lines
649
$line[0] = get_symbol($line[0]);
652
my $mess = join(" ", $save, $$messages{not_an_element},
653
$$messages{at_line}, $iline, $/);
657
##my $mess = join(" ", $save, "is not an element symbol", $/);
658
##$keys -> warn_or_die($mess, $die);
661
## ----- evaluate the coordinates safely
662
## see `perldoc Safe' for details
664
## this is a much weaker form of security, but it works with
666
my @xyz = map {number($_, $die, $keys)} ($line[1], $line[2], $line[3]);
668
## my $message = $$messages{'tainted_atoms'} . $/;
669
## $message .= $$messages{'file'} . ": $file, ";
670
## $message .= $$messages{'line'} . ": $.$/";
671
## my $cpt = new Safe;
672
## ##$cpt->share('@xyz');
673
## $ {$cpt->varglob('x')} = $cpt->reval($xyz[0]);
674
## $ {$cpt->varglob('y')} = $cpt->reval($xyz[1]);
675
## $ {$cpt->varglob('z')} = $cpt->reval($xyz[2]);
676
## @xyz = ($ {$cpt->varglob('x')}, $ {$cpt->varglob('y')},
677
## $ {$cpt->varglob('z')} );
678
## test_safe_return($message, @xyz);
679
# @xyz = map {sprintf "%11.7f", $_} @xyz;
680
## ----- end of safe evaluation
682
if (defined $line[4]) {
683
## if the fourth column matches a C float, assume it is an
684
## occupancy and that the tag should be ""
685
if ($line[4] =~ /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/) {
692
## test fifth column to see if it is an occupancy
693
(defined $line[5]) || ($line[5] = 1);
694
if ($line[5] =~ /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/) {
695
$line[5] = number($line[5], $die, $keys);
696
($line[5] =~ /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/)
701
$keys -> make('sites'=>
702
$line[0], $xyz[0], $xyz[1], $xyz[2], $line[4], $line[5]);
704
my $this = lc($line[0]);
707
## unknown keyword, store it anyway
709
my $value = shift @line;
710
$keys -> make($word=>$value);
711
## my $message = "\"$word\": " . $$messages{'unknown_keyword'} .
712
## ", " . $$messages{'line'} . "$..$/";
725
This method reads from a CIF file and loads the data into keywords
731
## this is ungainly due to the mismatch in nomenclature and verbosity
732
## of my stuff and the STAR stuff.
734
my ($keys, $file, $which) = @_;
735
my @datablocks = STAR::Parser->parse($file);
736
my $datablock = $datablocks[$which];
737
## if STAR::Checker->check(-datablock=>$datablocks[0]);
741
## titles: consider various common title-like entries, strip white
742
## space characters and stuff them into the title attribute
743
@item = $datablock->get_item_data(-item=>"_chemical_name_mineral");
744
$item[0] ||= ""; $item[0] =~ s/
745
$keys->make(title=>$item[0]) unless ($item[0] =~ /^\s*$/);
746
@item = $datablock->get_item_data(-item=>"_chemical_name_systematic");
747
$item[0] ||= ""; $item[0] =~ s/
748
$keys->make(title=>$item[0]) unless ($item[0] =~ /^\s*$/);
749
@item = $datablock->get_item_data(-item=>"_chemical_formula_structural");
750
$item[0] ||= ""; $item[0] =~ s/
751
$keys->make(title=>$item[0]) unless ($item[0] =~ /^\s*$/);
752
@item = $datablock->get_item_data(-item=>"_chemical_formula_sum");
753
$item[0] ||= ""; $item[0] =~ s/
754
$keys->make(title=>$item[0]) unless ($item[0] =~ /^\s*$/);
755
@item = $datablock->get_item_data(-item=>"_publ_author_name");
756
$item[0] ||= ""; $item[0] =~ s/
757
$keys->make(title=>$item[0]) unless ($item[0] =~ /^\s*$/);
758
@item = $datablock->get_item_data(-item=>"_citation_journal_abbrev");
759
$item[0] ||= ""; $item[0] =~ s/
760
$keys->make(title=>$item[0]) unless ($item[0] =~ /^\s*$/);
761
@item = $datablock->get_item_data(-item=>"_publ_section_title");
762
$item[0] ||= ""; $item[0] =~ s/
763
$keys->make(title=>$item[0]) unless ($item[0] =~ /^\s*$/);
765
## space group: try the number then the NM symbol and canonicalize it
766
@item = $datablock->get_item_data(-item=>"_symmetry_Int_Tables_number");
767
my @sg = Xray::Xtal::Cell::canonicalize_symbol($item[0]);
769
@item = $datablock->get_item_data(-item=>"_symmetry_space_group_name_H-M");
770
@sg = Xray::Xtal::Cell::canonicalize_symbol($item[0]);
772
$keys->make(space=>$sg[0]) if $sg[0];
775
my $min = 100000; # use lattice constants to compute default for Rmax
776
foreach my $k (qw(a b c)) {
777
@item = $datablock->get_item_data(-item=>"_cell_length_$k");
778
(my $this = $item[0]) =~ s/\(\d+\)//;
780
$keys->make($k=>$this);
781
$min = 1.1*$this if ($min > 1.1*$this);
783
$min = 7 if ($min > 11);
784
$keys->make(rmax=>$min);
785
foreach my $k (qw(alpha beta gamma)) {
786
@item = $datablock->get_item_data(-item=>"_cell_angle_$k");
787
(my $this = $item[0]) =~ s/\(\d+\)//;
789
$keys->make($k=>$this);
792
## load up and clean up the atom positions
793
my @tag = $datablock->get_item_data(-item=>"_atom_site_label");
794
my @el = $datablock->get_item_data(-item=>"_atom_site_type_symbol");
795
my @x = $datablock->get_item_data(-item=>"_atom_site_fract_x");
796
my @y = $datablock->get_item_data(-item=>"_atom_site_fract_y");
797
my @z = $datablock->get_item_data(-item=>"_atom_site_fract_z");
798
my @occ = $datablock->get_item_data(-item=>"_atom_site_occupancy");
799
foreach my $i (0 .. $#tag) {
800
my $ee = $el[$i] || $tag[$i];
802
(my $xx = $x[$i]) =~ s/\(\d+\)//; # remove parenthesized error bars
803
(my $yy = $y[$i]) =~ s/\(\d+\)//;
804
(my $zz = $z[$i]) =~ s/\(\d+\)//;
805
(my $oo = $occ[$i]||1) =~ s/\(\d+\)//;
806
##print "$ee, $xx, $yy, $zz, $tag[$i], $oo\n";
807
$keys -> make('sites'=> $ee, $xx, $yy, $zz, $tag[$i], $oo);
811
## Tags in the cif file seem to mostly be the two letter element
812
## symbol concatinated with a number (possibly concatinated with a
813
## non-alphanumeric, such as a plus sign or parens). Also "Wat" and "OH"
814
## indicate oxygen sites (with one or two associated hydrogens) and D
818
($elem =~ /Wat/) and return "O";
819
($elem =~ /OH/) and return "O";
820
($elem =~ /^D$/) and return "H";
821
## snip off the last character until an element symbol is found
823
return $elem if ($elem =~ /^$elem_match$/o);
830
sub test_safe_return {
831
my $message = shift(@_);
833
foreach my $item (@list) {
834
(defined $item) || die $message . $/;
839
=head2 C<verify_keywords>
841
This method is used to perform some sanity checks on the keyword
842
values after reading in the input data and populating the sites and the
843
cell. It is a good idea to call this before building a cluster or
844
starting some other calculation.
846
$keywords -> verify_keywords($cell, \@sites, 0, 0);
848
The first two arguments are a Cell object and a reference to a
849
list of Site objects. The third argument is the same as the second
850
argument to the C<parse_input> method. The last argument should be
851
non-zero only for a calculation that does not need to know which is
852
the central atom (e.g. a powder diffraction simulation).
856
sub verify_keywords {
858
my ($cell,$r_sites,$die,$no_core) = @_;
860
$keys ->{dopant_core} = "";
861
$keys -> set_rmax($cell, $def, $die);
863
my $checkit = $keys -> set_core($r_sites, $die);
864
return 1 if ($checkit eq '-1');
865
#$keys -> check_core($r_sites, $die);
866
$keys -> set_edge($cell,$die);
868
foreach my $word (qw/argon krypton nitrogen/) {
869
if (($keys->{$word} < 0) || ($keys->{$word} > 1)) {
870
$keys->make($word=>'0');
871
my $message = "\"$word\": $$messages{'one_gas'}$/";
872
$keys -> warn_or_die($message, $die);
875
if (($keys->{nitrogen} + $keys->{argon} + $keys->{krypton}) > 1) {
876
my $message = $$messages{'all_gases'} . $/;
877
$keys->make('nitrogen'=> '0');
878
$keys->make('argon'=> '0');
879
$keys->make('krypton'=> '0');
880
$keys -> warn_or_die($message, $die);
883
if ($keys->{dopant_core}) {
884
my $z = &get_Z($keys->{dopant_core});
886
my $message = $$::messages{dopant_core} . $/;
887
$keys -> warn_or_die($message, $die);
892
# perhaps take as arguments refs to subroutines for further tests
894
## the next several subroutines are used internally, thus not pod
897
## reference to a cell and default rmax value (defaults to 7)
900
my ($cell, $def, $die) = @_;
901
my $message = $$messages{'rmax'} . $/;
903
my $rm = number($keys->{'rmax'}, $die, $keys);
904
if (($rm != 0) && ($rm < EPSILON)) {
906
$keys -> warn_or_die($message, $die);
908
($rm > EPSILON) || do {
909
## rmax is the larger of (7 and 1.1 times the shortest axis)
910
my @list = sort ($cell->attributes('A','B','C'));
911
@list = sort ($list[0], $def/1.1);
912
$keys->{"rmax"} = 1.1*$list[1];
914
return $keys->{"rmax"};
917
## reference to an array of sites
922
#unless ($keys->{"core"}) {
923
if ($#{$sites} == 0) {
924
my ($t) = $$sites[0] -> attributes('Tag');
925
$keys->make("core", $t);
926
return $keys->{"core"};
927
} elsif (($#{$sites} != 0) and ($keys->{"core"} =~ /^\s*$/)) {
928
$keys -> warn_or_die($$messages{'no_core'} . $/, $die);
929
my ($t) = $$sites[0] -> attributes('Tag');
930
$keys->make("core", $t);
931
return $keys->{"core"};
932
} elsif ($keys -> check_core($sites, $die)) {
933
return $keys->{"core"};
935
$keys -> warn_or_die($$messages{'no_core'} . $/, $die);
937
#return $keys->{"core"};
943
my ($sites, $die) = @_;
944
my $c = lc($keys->{'core'});
946
foreach my $site (@$sites) {
947
my ($t) = $site -> attributes('Tag');
949
$found = 1, last if ($t eq $c);
953
$keys -> warn_or_die("\"" . $keys->{'core'} . "\" " .
954
$$messages{'unknown_core'} . $/,
956
(not $found) and ($die == 0) and die "\n";
957
#$::help{'check_core'};
963
my ($cell,$die) = @_;
964
if ($cell->{Contents}) {
965
unless ($keys->{'edge'}) { # default based on Z number
966
my ($central, $xcenter, $ycenter, $zcenter) =
967
$cell -> central($keys->{"core"});
968
my $z = &get_Z($central);
969
$keys->{'edge'} = 'K';
971
$keys->{'edge'} = 'L3';
975
($keys->{"edge"} =~ /k|l[123]|m[1-5]|n[1-7]|o[1-7]|p[1-3]/i) || do {
976
my $message = $$messages{'bad_edge'} . $/;
977
($keys->{"edge"} = 'k');
978
$keys -> warn_or_die($message, $die);
982
return $keys->{'edge'};
989
return $keys->{'atp'}{$atp};
994
if ($keys->{'atp'}{$atp}) {
995
$keys->{'atp'}{$atp} = 0;
997
$keys->{'atp'}{$atp} = 1;
1001
=head2 C<warn_or_die>
1003
This is a subroutine rather than a method, so that is can be used even
1004
without an active keywords object. It is a wrapper for generating
1005
error and warning methods that will work in command-line, GUI, and CGI
1008
$keywords -> warn_or_die("This is a warning message", $die);
1010
The first argument is the actual message. The C<die> argument is 0,
1011
1, or 2 as described for C<parse_input>. The final argument is a
1012
keywords object and is only used in the CGI environment. In that
1013
case, here are no convenient output or error channels, so messages are
1014
stored in the keywords hash as a single string under the key
1021
my ($message, $die) = @_;
1022
if ($die == 1) { # Tk
1023
##my $temp = MainWindow->new;
1024
##$temp -> withdraw;
1025
$::top -> messageBox(-icon => 'error',
1026
-message => $message,
1027
-title => 'Atoms: Error',
1031
} elsif ($die == 2) { # CGI
1032
$keys -> {www_warn} .= $message;
1033
} elsif ($die == 0) { # CLI
1035
$keys -> {cli_warn} .= $message if $keys;
1036
} elsif ($die == 4) { # ignore
1039
die "Programmer error! Invalid run level $die in warn_or_die!\n";
1043
## =head2 C<available_languages>
1045
## This subroutine takes no arguments and returns a list of the languages
1046
## for which translations of the language data exist.
1050
sub available_languages {
1051
return @available_languages;
1057
This subroutine is the output engine. It parses the appropriate atp
1058
file, interprets its content, and generates output data. This one
1059
subroutine handles the creation of F<feff.inp> and the other sorts of
1060
lists generated by Atoms. It also generates reports for the
1061
Absorption notecard in TkAtoms and for the lists describing the DAFS
1062
calculation as a function of energy and the powder calculation as a
1063
function of angle. Essentially all output from all programs that come
1064
with the Atoms package is generated here.
1066
my ($default_name, $is_feff) =
1067
&parse_atp($atp, $cell, $keywords, \@cluster, \@neutral, \$contents);
1068
open (INP, ">".$default_name) or die $!;
1069
print INP $contents;
1078
A string specifying the atp file to read. The actual atp file either
1079
lives in the F<atp/> subdirectory of the Atoms installation or in the
1080
F<.atoms/> directory of the individual user.
1092
A reference to the list generated by the C<build_cluster> subroutine,
1093
to the list containing a calculation such as the one in DAFS or
1094
Powder, or 'whatever the thing is that I did to make atoms.inp output
1099
A reference to a list. This is just a place holder and is not
1100
currently used for anything.
1104
A scalar reference. This is assumed to be empty on entry, and is
1105
filled with contents intended for the output file on exit.
1109
Typically the C<contents> variables is the things written to the
1110
output channel. The two return values are the default name for the
1111
output file and a flag indicating whether the output is intended to
1112
run FEFF. Both of these pieces of information are read from the
1113
C<meta> line of the atp file.
1117
=head1 EXPORTED SUBROUTINES
1119
=head2 C<absorption>
1121
This subroutine calculates the total absorption (or one over the
1122
e-fold absorption length), the delta_mu across the edge (or one over
1123
the unit edge step length) and the specific gravity of the crystal.
1125
($total_absorption, $delta_mu, $specific_gravity) =
1126
&absorption(\$cell, $central, $edge);
1128
The input arguments are a populated cell, the tag of the central atom,
1129
and the alphanumeric (i.e. K, L1, etc) symbol of the absorption edge.
1131
Note that the values returned depend on the data resource used. See
1132
L<Xray::Absorption>.
1138
## this reproduces well the results from the Fortran atoms.
1140
my ($cell, $central, $edge) = @_;
1141
$central = lc($central);
1142
my $specified = ($edge =~ /\b\d+\b/);
1143
my $energy = ($edge =~ /\b\d+\b/) ?
1144
$edge : Xray::Absorption -> get_energy($central, $edge);
1145
my ($contents) = $cell -> attributes("Contents");
1146
my ($bravais) = $cell -> attributes("Bravais");
1147
my $brav = ($#{$bravais}+4) / 3;
1148
my ($volume, $occ) = $cell -> attributes("Volume", "Occupancy");
1149
my ($mass, $xsec, $delta_mu, $den, $conv) = (0,0,0,0,0);
1150
my %cache = (); # memoize and call cross_section less often
1151
foreach my $atom (@{$contents}) {
1152
my ($element, $this_occ) =
1153
$ {$$atom[3]} -> attributes("Element", "Occupancy");
1154
# print join(" ", $element, $this_occ, $/);
1155
$element = lc($element);
1156
my $factor = $this_occ; # $occ ? $this_occ : 1; # consider site occupancy??
1157
my $weight = Xray::Absorption -> get_atomic_weight($element);
1158
$mass += $weight*$factor;
1159
$cache{lc($element)} ||=
1160
scalar Xray::Absorption -> cross_section($element, $energy+50);
1161
$xsec += $cache{lc($element)} * $factor;
1162
if (($central eq $element) and not $specified) {
1165
( $cache{lc($central)} -
1166
scalar Xray::Absorption -> cross_section($central, $energy-50) );
1169
$mass *= 1.66053/$volume; ## atomic mass unit = 1.66053e-24 gram
1171
$delta_mu /= $volume;
1172
return ($xsec, $delta_mu, $mass);
1178
This subroutine calculates the total absorption (or one over the
1179
e-fold absorption length) and the specific gravity of the crystal at a
1180
specified energy (rather than around the edge energy of an element).
1184
($total_absorption, $specific_gravity) =
1185
&xsec($cell, $central, $energy);
1189
(\@total_absorption, $specific_gravity) =
1190
&xsec($cell, $central, \@energies);
1192
The input arguments are a populated cell and the tag of the central
1193
atom. The third argument is either the desired energy in eV or a
1194
reference to a list of energies in eV.
1196
Note that the values returned depend on the data resource used. See
1197
L<Xray::Absorption>.
1203
my ($cell, $central, $energy) = @_;
1204
$central = lc($central);
1205
my ($contents) = $cell -> attributes("Contents");
1206
my ($bravais) = $cell -> attributes("Bravais");
1207
my $brav = ($#{$bravais}+4) / 3;
1208
my ($volume, $occ) = $cell -> attributes("Volume", "Occupancy");
1209
my ($xsec, $delta_mu, $den, $conv) = (0,0,0,0);
1211
my %cache = (); # memoize and call cross_section less often
1212
foreach my $atom (@{$contents}) {
1213
my ($element, $this_occ) =
1214
$ {$$atom[3]} -> attributes("Element", "Occupancy");
1215
$element = lc($element);
1216
my $factor = $this_occ; # $occ ? $this_occ : 1; # consider site occupancy??
1218
exists $cache{$element} or
1219
@{$cache{$element}} =
1220
Xray::Absorption -> cross_section($element, $energy);
1221
foreach (0 .. $#{$energy}) {
1222
$xsec[$_] += $ {$cache{$element}}[$_] * $factor;
1225
$cache{$element} ||=
1226
scalar Xray::Absorption -> cross_section($element, $energy);
1227
$xsec += $cache{$element} * $factor;
1231
return wantarray ? @xsec : $xsec;
1236
my ($contents, $bravais) = $cell -> attributes("Contents", "Bravais");
1237
my $brav = ($#{$bravais}+4) / 3;
1238
my ($volume, $occ) = $cell -> attributes("Volume", "Occupancy");
1239
my ($mass, $den, $conv) = (0,0,0);
1240
foreach my $atom (@{$contents}) {
1241
my ($element, $this_occ) =
1242
$ {$$atom[3]} -> attributes("Element", "Occupancy");
1243
$element = lc($element);
1244
my $factor = $this_occ; # $occ ? $this_occ : 1; # consider site occupancy??
1245
my $weight = Xray::Absorption -> get_atomic_weight($element);
1246
$mass += $weight*$factor;
1248
$mass *= 1.66053/$volume; ## atomic mass unit = 1.66053e-24 gram
1253
sub mcmaster_pre_edge {
1254
my ($central, $edge) = @_;
1256
my $emin = Xray::Absorption -> get_energy($central, $edge) - 10;
1257
## find the pre-edge line
1258
my %next_e = ("k"=>"l1", "l1"=>"l2", "l2"=>"l3", "l3"=>"m");
1260
if (exists $next_e{$edge}) {
1261
$ebelow = Xray::Absorption -> get_energy($central, $next_e{$edge}) + 10;
1262
$ebelow = (($emin - $ebelow) > 100) ? $emin - 100 : $ebelow;
1264
$ebelow = $emin - 100;
1266
my $delta = ($emin - $ebelow)/10;;
1267
my @i=(0..9); # load the pre edge energies/sigmas
1268
my @energy = map {$ebelow + $delta*$_} @i;
1269
my @sigma = Xray::Absorption -> cross_section($central, \@energy);
1271
my $pre_edge = Statistics::Descriptive::Full->new();
1272
$pre_edge -> add_data(@sigma);
1273
my ($bpre, $slope) = $pre_edge -> least_squares_fit(@energy);
1274
$bpre ||= 0; $slope ||= 0;
1275
return ($bpre, $slope);
1280
This is called C<mcmaster> for historical reasons. It calculates the
1281
normalization correcion for a given central atom.
1283
$sigma_mm = &mcmaster($central, $edge);
1285
It takes the central atoms tag and the alphanumeric edge symbol as
1286
arguments and returns the normalization correction in units of
1289
Note that the values returned depend on the data resource used. See
1290
L<Xray::Absorption>.
1295
## $span in the fortran version was 500eV and the regression was
1296
## performed with a square term.
1297
## Statistics::Descriptive::least_squares_fit only fits a line, so I
1298
## drew the #span back to 300 volts. This gives the "canonical"
1299
## 0.00052 for copper.
1301
my ($central, $edge) = @_;
1303
my $emin = Xray::Absorption -> get_energy($central, $edge) + 10 ;
1304
my ($emax, $e_to_k, $span, $npost) = (0, 0.2624682917, 300, 20);
1306
## get energy range for above edge fit
1307
my %next_e = ("l1"=>"k", "l2"=>"l1", "l3"=>"l2", "m"=>"l3");
1308
if (exists $next_e{$edge}) {
1309
$emax = Xray::Absorption -> get_energy($central, $next_e{$edge}) - 10;
1310
$emax = (($emax - $emin) > $span) ? $emin + $span : $emax;
1312
$emax = $emin + $span;
1315
## need to show some care with the Chantler data
1316
(Xray::Absorption -> current_resource() =~ /chantler/i) and
1318
($emin >= $emax) and ($emin = $emax - 20); # whatever!
1320
my ($bpre, $slope) = mcmaster_pre_edge($central, $edge);
1321
my $delta = ($emax - $emin)/$npost;
1322
my @i=(0..$npost-1); # load the post edge energies and sigmas
1323
my @energy = map {$emin + $delta*$_} @i;
1324
## and some more care...
1325
(Xray::Absorption -> current_resource() =~ /chantler/i) and do {
1326
shift @energy; shift @energy;
1328
return 0 if ($bpre <= 0);
1329
my @sigma = Xray::Absorption -> cross_section($central, \@energy);
1330
@sigma = map {$sigma[$_] - ($bpre+$energy[$_]*$slope)} (0 .. $#energy);
1331
##map {printf " %9.3f %9.3f\n", $energy[$_], $sigma[$_]} (0 .. $#energy);
1332
@energy = map {$e_to_k * ($_-$emin)} @energy; # convert to k
1333
my $any_neg = grep {$_ <= 0} @sigma;
1334
return 0 if $any_neg;
1335
@sigma = map {log($_)} @sigma; # take logs of xsecs
1337
my $stat = Statistics::Descriptive::Full->new(); # fit the post edge
1338
$stat -> add_data(@sigma);
1339
my @a = $stat -> least_squares_fit(@energy);
1340
return ($a[1] < 0) ? -$a[1]/2 : 0;
1346
This calculates the correcion due to the I0 fill gases in a
1347
fluorescence experiment.
1349
$sigma_i0 = &i_zero($central, $edge, $nitrogen, $argon, $krypton);
1351
It takes the central atoms tag, the alphanumeric edge symbol, and the
1352
volume percentages of the three gases as arguments. It assumes that
1353
any remaining volume is filled with helium and it correctly accounts
1354
for the fact that nitrogen is a diatom. It returns the I0 correction
1355
in units of Angstrom squared.
1357
Note that the values returned depend on the data resource used. See
1358
L<Xray::Absorption>.
1363
my ($central, $edge, $nitrogen, $argon, $krypton) = @_;
1365
## convert from pressure percentages to number of absorbers.
1366
## nitrogen is diatomic
1367
my $helium = 1 - $nitrogen - $argon - $krypton;
1368
my $norm = $helium + 2*$nitrogen + $argon + $krypton;
1369
$nitrogen = 2*$nitrogen / $norm;
1370
$argon = $argon / $norm;
1371
$krypton = $krypton / $norm;
1373
my $emin = Xray::Absorption -> get_energy($central, $edge) ;
1374
my ($emax, $e_to_k, $span, $npost) = (0, 0.2624682917, 500, 20);
1375
## careful not to run a gas edge
1376
my ($el, $ed, $en) = Xray::Absorption ->
1377
next_energy($central, $edge, "ar", "n", "kr");
1378
if (not defined $en) {
1379
$emax = $emin + $span;
1381
$emax = ( ($en - 10) < ($emin + $span) ) ? ($en - 10) : ($emin + $span);
1383
## need to show some care with the Chantler data
1384
(Xray::Absorption -> current_resource() =~ /chantler/i) and
1386
($emin >= $emax) and ($emin = $emax - 20); # whatever!
1388
my @i=(0..$npost-1); # load the post edge energies and sigmas
1389
my $delta = ($emax - $emin)/$npost;
1390
my @energy = map {$emin + $delta*$_} @i;
1391
my @s_n = Xray::Absorption -> cross_section("n", \@energy);
1392
my @s_a = Xray::Absorption -> cross_section("ar", \@energy);
1393
my @s_k = Xray::Absorption -> cross_section("kr", \@energy);
1395
{$nitrogen*$s_n[$_] + $argon*$s_a[$_] + $krypton*$s_k[$_]} (0 .. $#energy);
1396
@energy = map {$e_to_k * ($_-$emin)} @energy; # convert to k
1397
@sigma = map {log($_)} @sigma; # take logs of xsecs
1399
my $stat = Statistics::Descriptive::Full->new(); # fit the post edge
1400
$stat -> add_data(@sigma);
1401
my @a = $stat -> least_squares_fit(@energy);
1407
This calculates the correcion due to self-absorption fluorescence
1408
experiment. It assumes that the sample is infinately thick and that
1409
the entry and exit angles of the photons are the same.
1411
$sigma_i0 = &self($central, $edge, $cell);
1413
It takes the central atoms tag, the alphanumeric edge symbol, and a
1414
populated cell. It returns a list whose zeroth element is the
1415
multiplicative amplitude correction and whose first element is the a
1416
correction in units of Angstrom squared.
1418
Note that the values returned depend on the data resource used. See
1419
L<Xray::Absorption>.
1424
my ($central, $edge, $cell) = @_;
1426
my ($contents) = $cell -> attributes("Contents");
1427
my ($occ) = $cell -> attributes("Occupancy");
1428
foreach my $atom (@{$contents}) {
1429
push @list, $ {$$atom[3]} -> attributes("Element");
1432
my $emin = Xray::Absorption -> get_energy($central, $edge) ;
1433
my ($emax, $e_to_k, $span, $npost) = (0, 0.2624682917, 500, 20);
1434
my ($el, $ed, $en) = Xray::Absorption -> next_energy($central, $edge, @list);
1435
if (not defined $en) {
1436
$emax = $emin + $span;
1438
$emax = ( ($en - 10) < ($emin + $span) ) ? ($en - 10) : ($emin + $span);
1440
## need to show some care with the Chantler data
1441
(Xray::Absorption -> current_resource() =~ /chantler/i) and
1443
($emin >= $emax) and ($emin = $emax - 20); # whatever!
1445
## calculate total absorption at the fluorescence energy and 10
1446
## volts below the edge
1448
my $fline = substr($edge, 0, 1) . "alpha";
1449
my $e_fluor = Xray::Absorption -> get_energy($central, $fline);
1450
#my $e_below = $emin - 10;
1451
foreach my $atom (@{$contents}) {
1452
my ($element, $this_occ) = $ {$$atom[3]} -> attributes("Element", "Occupancy");
1453
my $factor = $this_occ ; # $occ ? $this_occ : 1; # consider site occupancy??
1454
$xmuf += scalar Xray::Absorption -> cross_section($element, $e_fluor) * $factor;
1455
#$xmub += cross_section($element, $e_below);
1458
## load the post edge energies and sigmas
1461
my @i=(0..$npost-1);
1462
my $delta = ($emax - $emin)/$npost;
1463
my @energy = map {$emin + $delta*$_} @i;
1466
foreach my $j (@i) {
1469
my %cache = (); # memoize and call cross_section less often
1470
foreach my $atom (@{$contents}) {
1471
my ($element, $this_occ) = $ {$$atom[3]} -> attributes("Element", "Occupancy");
1472
if (lc($element) eq lc($central)) {
1473
my $factor = $this_occ; # $occ ? $this_occ : 1;
1474
$cache{lc($element)} ||=
1475
Xray::Absorption -> cross_section($element, $energy[$j]) * $factor;
1476
$xmu_core += $cache{lc($element)};
1478
my $factor = $this_occ;
1479
$cache{lc($element)} ||=
1480
Xray::Absorption -> cross_section($element, $energy[$j]) * $factor;
1481
$xmu += $cache{lc($element)};
1484
$sigma[$j] = ($xmuf+$xmu+$xmu_core)/($xmuf+$xmu);
1487
@energy = map {$e_to_k * ($_-$emin)} @energy; # convert to k
1488
@sigma = map {log($_)} @sigma;
1490
my $stat = Statistics::Descriptive::Full->new(); # fit the post edge
1491
$stat -> add_data(@sigma);
1492
my @a = $stat -> least_squares_fit(@energy);
1493
return (exp($a[0]), -$a[1]/2);
1499
## =item &bravais_string
1501
## Returns the peculiar array from &bravais as a pretty-printed string.
1502
## This is not method of the Cell class, it is just a normal function.
1503
## This demonstrates the differnt return values of bravais and
1506
## @list = Xray::Xtal::Cell::bravais("F m 3 m", 0);
1507
## print Xray::Xtal::Cell::bravais_string(@list), "\n",
1508
## join(", ", @list), "\n":
1510
## |-> (0, 0, 0), (0, 1/2, 1/2), (1/2, 0, 1/2), (1/2, 1/2, 0)
1511
## (0, 0, 0, 0, 0.5, 0.5, 0.5, 0, 0.5, 0.5, 0.5, 0)
1517
sub bravais_string {
1518
my ($bravais, $gnxas) = @_;
1519
return "" unless ($bravais);
1520
my @temp = @$bravais;
1521
length(join("", @temp)) || return ($gnxas) ? "1$/0,0,0" : "(0, 0, 0)";
1523
if (abs($_ - 0.5) < EPSILON) {
1525
} elsif (abs($_ - 2/3) < EPSILON) {
1527
} elsif (abs($_ - 1/3) < EPSILON) {
1535
($#temp == 2) and $string = "2$/";
1536
($#temp == 5) and $string = "3$/";
1537
($#temp == 8) and $string = "4$/";
1538
$string .= "0.0000,0.0000,0.0000";
1540
$string .= sprintf("$/%6.4f,%6.4f,%6.4f",
1541
eval(shift @temp), eval(shift @temp), eval(shift @temp));
1544
$string = "(0, 0, 0), (";
1545
$string .= join(", ", $temp[0], $temp[1], $temp[2]);
1547
($string .= "), (" . join(", ", $temp[3], $temp[4], $temp[5]));
1549
($string .= "), (" . join(", ", $temp[6], $temp[7], $temp[8]));
1555
=head2 C<build_cluster>
1557
This exported routine builds a spherical cluster from a populated unit
1558
cell by constructing a rhomboid of an intger number of complete cells
1559
which fully encloses the sphere. The spherical cluster is probably
1562
build_cluster($cell, $keywords, \@cluster, \@neutral);
1564
The first two arguments are cell and keywords objects. The cell
1565
should already have been populated by the C<populate> method. The
1566
thrid argument is a reference to an array that will contain the
1567
cluster. The fourth argument is an unused placeholder.
1569
The C<@cluster> array is an anonymous array or arrays. Each entry in
1570
C<@cluster> is an array containing x, y, z, r, r squared, x, y, z, the
1571
formula for x, the formula for y , and the formula for z. The first
1572
three arguments are the coordinates to full precision. The fourth
1573
through seventh arguments are of limited precision and are used
1574
internally for sorting the cluster.
1579
# cell, rmax, position of central atom
1580
# return cluster and neutral cluster
1583
my ($cell, $keys, $r_cluster, $r_neutral) = @_;
1584
my $rmax = $keys->{'rmax'};
1587
my ($central, $xcenter, $ycenter, $zcenter) =
1588
$cell -> central($keys->{"core"});
1589
my $setting = $cell->{Setting};
1590
my $crystal_class = $cell -> crystal_class();
1591
my $do_tetr = ($crystal_class eq "tetragonal" ) && ($setting);
1592
##print join(" ", $central, $xcenter, $ycenter, $zcenter), $/;
1593
my ($aa, $bb, $cc) = $cell -> attributes("A", "B", "C");
1594
my $xup = ceil($rmax/$aa - 1 + $xcenter);
1595
my $xdn = ceil($rmax/$aa - $xcenter);
1596
my $yup = ceil($rmax/$bb - 1 + $ycenter);
1597
my $ydn = ceil($rmax/$bb - $ycenter);
1598
my $zup = ceil($rmax/$cc - 1 + $zcenter);
1599
my $zdn = ceil($rmax/$cc - $zcenter);
1600
##print join(" ", "up,dn", $xup, $xdn, $yup, $ydn, $zup, $zdn), $/;
1601
#my $num_z = int($rmax/$cc) + 1; # |
1602
my $rmax_squared = $rmax**2; # (sprintf "%9.5f", $rmax**2);
1603
my ($contents) = $cell -> attributes("Contents");
1604
foreach my $nz (-$zdn .. $zup) {
1605
foreach my $ny (-$ydn .. $yup) {
1606
foreach my $nx (-$xdn .. $xup) {
1607
foreach my $pos (@{$contents}) {
1608
my ($x, $y, $z) = ($$pos[0]+$nx, $$pos[1]+$ny, $$pos[2]+$nz);
1609
($x, $y, $z) = ($x-$xcenter, $y-$ycenter, $z-$zcenter);
1610
($x, $y, $z) = $cell -> metric($x, $y, $z);
1611
($do_tetr) and ($x, $y) = (($x+$y)/sqrt(2), ($x-$y)/sqrt(2));
1612
#printf "in: %25s %25s %25s %3d %3d %3d\n", @$pos[4..6],$nx,$ny,$nz;
1613
my ($fx, $fy, $fz) = &rectify_formula(@$pos[4..6], $nx, $ny, $nz);
1614
#printf "out: %25s %25s %25s\n\n", $fx, $fy, $fz;
1615
my $r_squared = (sprintf "%9.5f", $x**2 + $y**2 + $z**2);
1616
my $this_site = [$x, $y, $z, $$pos[3],
1617
$r_squared, # cache the
1618
(sprintf "%11.7f", $x), # stuff needed
1619
(sprintf "%11.7f", $y), # for sorting
1620
(sprintf "%11.7f", $z),
1622
($r_squared < $rmax_squared) && # keep the ones within rmax
1623
(push @$r_cluster, $this_site);
1624
## (push @$r_neutral, $this_site);
1630
## =============================== sort the cluster (& neutral clus.)
1631
foreach my $r_list ($r_cluster) { ##, $r_neutral) {
1633
($a->[4] cmp $b->[4]) # sort by distance squared or ...
1635
($ {$a->[3]}->{Element} cmp $ {$b->[3]}->{Element})
1637
($a->[7] cmp $b->[7]) # by z value or ...
1639
($a->[6] cmp $b->[6]) # by y value or ...
1641
($a->[5] cmp $b->[5]); # by x value
1643
## ($ {$b->[3]}->{Host} <=> $ {$a->[3]}->{Host}); # hosts before dopants
1647
## final adjustment to the formulas, store the formulas for the
1649
$keys -> {cformulas} =
1650
[$$r_cluster[0][8], $$r_cluster[0][9], $$r_cluster[0][10]];
1651
## ## ... subtract the central atom coordinates from each site ...
1652
## foreach my $site (reverse(@$r_cluster)) {
1653
## (@$site[8..10]) =
1654
## ($$site[8] . " - Xc", $$site[9] . " - Yc", $$site[10] . " - Zc");
1656
## ## ... and set the central atom to an empty string
1657
## ($$r_cluster[0][8], $$r_cluster[0][9], $$r_cluster[0][10]) = ("", "", "");
1659
## if this is a tetragonal crystal in the C or F setting , rotate
1660
## all the coordinates back to the original setting
1662
my ($a, $b) = ($cell->{A}, $cell->{B});
1663
$cell -> make(A=>$a*sqrt(2), B=>$b*sqrt(2));
1668
sub rectify_formula {
1669
my ($fx, $fy, $fz, $nnx, $nny, $nnz, $r_cformula) = @_;
1670
my $debug_formulas = 1;
1671
##open DEBUG, '>debug.formulas' or die $!;
1672
#warn "input: ", join(" | ",$fx, $fy, $fz, $nnx, $nny, $nnz, $/);
1673
## apply bravais translation
1675
if ($_ < 0) { $_ = ' + '.$_ }
1676
elsif ($_ == 0) { $_ = '' }
1677
else { $_ = ' + '.$_ };
1678
} ($nnx, $nny, $nnz);
1679
($fx, $fy, $fz) = ($fx . $nnx, $fy . $nny, $fz . $nnz);
1680
## convert fractions to decimal
1681
map { s|(\d)/(\d)|$1/$2|ge } ($fx, $fy, $fz);
1682
# collapse the numeric fields
1684
my @list = reverse(split(/\s*\+\s*/, $_));
1685
#print join('|', @list, $/);
1687
LIST: while (@list) {
1688
my $this = shift @list;
1690
if ($this =~ /^-?[xyz]/) { # whoops! put the variable back on
1691
unshift @list, $this;
1694
$sum = eval '$sum + $this'; # eval the constants
1696
if ($sum) { # prettify the number and unshift it
1697
$sum = sprintf("%10.6f", $sum);
1700
unshift @list, $sum;
1702
@list = reverse(@list);
1703
$_ = join("+", @list); # make a math expression
1704
s/\+-/-/g; # finally clean up the + and - signs
1709
#warn "output: ", join(" ", $fx, $fy, $fz, $/ x 2);
1710
return ($fx, $fy, $fz);
1713
=head2 C<rcfile_name>
1715
This takes no arguments and returns the name of the Atoms runtime
1716
configuration file belonging to the user. This does the ``right
1717
thing'' on the different platforms.
1722
return "???" if (&rcdirectory() eq '???');
1723
return File::Spec -> catfile(&rcdirectory(), "tkatoms.ini")
1724
if (($^O eq 'MSWin32') or ($^O eq 'cygwin'));
1725
return File::Spec -> catfile(&rcdirectory(), "atomsrc");
1730
} elsif ($^O eq 'os2') {
1732
} elsif ($^O eq 'MacOS') {
1734
} elsif (($^O eq 'MSWin32') or ($^O eq 'cygwin')) {
1735
return File::Spec->catfile($ENV{FEFFIT_DIR}, "horae") if defined($ENV{FEFFIT_DIR});
1739
eval '$home = $ENV{"HOME"} || $ENV{"LOGDIR"} || (getpwuid($<))[7];' or
1741
return $home . "/.horae/";
1746
open RC, $_[0] or die "could not open $_[0] as a configuration file\n";
1757
unless ($_ =~ /\$c_/) {s/\s*\#.*$//;};
1758
my @line = split(/\s*[ \t=]\s*/, $_);
1759
(defined $line[1]) and
1760
(($line[1] eq "''") or ($line[1] eq '""')) and $line[1] = '';
1762
($line[0] =~ /^\s*\$([a-zA-Z_]+)/) and do {
1764
$line[1] =~ s/[;\'\"]//g;
1765
if ($var =~ /^c_/) { # colors
1766
## my $v = substr($var, 2);
1767
## $colors{$v} = $line[1];
1768
## if (($colors{$v} =~ /^[0-9a-fA-F]{6}$/) or
1769
## ($colors{$v} =~ /^[0-9a-fA-F]{12}$/)) {
1770
## $colors{$v} = '#'.$line[1];
1772
} elsif ($var =~ /^f_/) { # fonts
1773
## my $v = substr($var, 2);
1774
## $fonts{$v} = join(" ", @line[1..$#line]);
1775
## $fonts{$v} =~ s/[;\'\"]//g;
1777
$meta{$var} = $line[1];
1781
($mode eq 'meta') and do {
1782
$meta{$line[0]} = $line[1];
1785
### don't need to parse colors and fonts here
1786
## ($mode eq 'colors') and do {
1787
## $colors{$line[0]} = $line[1];
1788
## (($colors{$line[0]} =~ /[0-9a-fA-F]{6}/) or
1789
## ($colors{$line[0]} =~ /[0-9a-fA-F]{12}/)) and
1790
## $colors{$line[0]} = '#'.$line[1];
1793
## ($mode eq 'fonts') and do {
1794
## $fonts{$line[0]} = join(" ", @line[1..$#line]);
1803
This takes no argmuents and returns the hash containing the values of
1804
variables (but not of fonts or colors) read from the rc file.
1812
## return ($always_write_feff, $atoms_language, $write_to_pwd,
1814
## $prefer_feff_eight, $absorption_tables, $dafs_default,
1816
## $plotting_hook, $default_filepath, $display_balloons,
1818
## $no_crystal_warnings, $one_frame, $convolve_dafs,
1820
## $never_ask_to_save, $ADB_location)
1825
This takes a text string and attempts to evaluate it as a number. It
1826
uses eval and so allows for simple math expressions, but it tries to
1827
be safe and not eval just any old math expression. You can use this
1828
if you want to evaluate numbers in the same manner as Atoms.
1831
printf "%7.5f\n", $x;
1836
## this returns a text string as a proper number
1839
(my $input = $num) =~ s/\s//g; # trim blanks
1840
my $die = $_[1] || 0;
1841
my $keys = $_[2] || Xray::Atoms->new();
1843
my $num_match = '([+-]?(\d+\.?\d*|\.\d+))';
1844
($num =~ /^\s*$/) && return 0; # null value
1846
if ($num =~ /^\s*$num_match\s*$/) { # floating point number
1847
($num = sprintf "%9.5f", $num) =~ s/^\s+//;
1849
}; # simple binary operation
1850
if ($num =~ /^\s*($num_match)\s*(\-|\+|\/|\*)\s*($num_match)\s*$/) {
1852
($num = sprintf "%9.5f", $num) =~ s/^\s+//;
1855
#(abs($num) < $Xray::Atoms::epsilon) && ($num = 0);
1856
## interpret angle as 80.23'45"
1857
my $message = join("", "The string \"", $num, "\" was found among the atom coordinates. Atoms doesn't know what to do with it. The parameter was set to zero. You might want to verify your input data.", $/);
1858
$keys -> warn_or_die($message, $die);
1863
# Autoload methods go after =cut, and are processed by the autosplit program.
1870
=head1 MORE INFORMATION
1872
There is more information available in the Atoms document. There you
1873
will find complete descriptions of atp files, calculations using the
1874
Xray::Absorption package, keywords in atoms input files and lots of
1880
Bruce Ravel <ravel@phys.washington.edu>
1881
Atoms URL: http://feff.phys.washington.edu/~ravel/software/atoms/
1887
## time-stamp-line-limit: 25
added
removed
Atoms/Atoms.pm
