Новосибирский институт органической химии им. Н.Н. Ворожцова СО РАН
Лаборатория изучения механизмов органических реакций
 |
 |
 |
 |
 |
 |
g2pri_basis
#!/usr/bin/perl -ws #use Data::Dump 'pp'; our ($atom,$no_del,$gen,$file,$h,$help); $| = 1; (my $script = $0) =~ s/.*\///; if ($h || $help) { print <<HELP; Генерирование базисов (включая auxiliary) для Природы с помощью orca или g09 Usage: $script [-atom=1,4..9] basis_name Зависимости: perl, orca или g09 Базисы серии L достаточны для решения многих задач, но они нигде, кроме как в Природе не используются. Между тем, изредка возникает потребность в других, более стандартных базисах. Например, чтобы воспроизвести расчет, сделанный в другой программе. Этот скрипт генерирует для Природы базисы, доступные как стандартные для программ orca или g09. Эти программы запускается для каждого из элементов с ключевыми словами ORCA - '! PBE basis_name PrintBasis NoIter' гауссиан - '# test BLYP/basis_name DenFit GFInput Guess=Only Pop=None' Опции: -gen=orca|g09 Программа-генератор базисов. Default -gen=orca (генерируемые ею дополнительные базисы оказались лучше для Природы). -atom=1..104 Номера элементов, для которых нужны базисы. Default -atom='1,4..9' -no_del Не удалять инпуты/аутпуты orca и g09 (для тестирования) -file=file_name Названия базисов берутся из файла (бьется по пробелам/строкам). HELP exit; } $atom ||= '1,4..9'; $gen ||= 'orca'; if ($gen !~ /^(g09|orca)$/) { die "-gen=g09|orca\n"; } if ($file && -f $file) { open F, '<', $file or die "Can't open $file: $!\n"; my @l = <F>; close F; foreach (@l) { push @ARGV, split; } } die "usage $script [-atom=1,4..9] basis_name\n" unless @ARGV; my $g09 = 'rung09'; die "No $g09\n" if $gen eq 'go9' && ! -x $g09; my $orca = '/usr/local/bin/orca'; die "No $orca\n" if $gen eq 'orca' && ! -x $orca; # Вгоняем таблицу Менделеева my @ATOM = qw (XX H He Li Be B C N O F Ne Na Mg Al Si P S Cl Ar K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe Cs Ba La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn Fr Ra Ac Th Pa U Np Pu Am Cm Bk Cf Es Fm Md No Lr Rf Db Sg Bh Hs Mt Ds Rg ); my %SPD = ('S'=>0,'P'=>1,'D'=>2,'F'=>3,'G'=>4,'H'=>5, 'I'=>6,'J'=>7,'K'=>8,'L'=>8,'M'=>10,'N'=>11); BASIS: foreach my $basis_name (@ARGV) { print "$basis_name: "; my $info = " Basis and auxiliary basis sets for PRIRODA generated by $script script from $gen program outputs\n"; $info .= sprintf(" %10s%-30s%s\n",'','basis','basis2') if $gen eq 'orca'; my $basis = "\$basis\ntype=general\ndefault=$basis_name\n\nset=$basis_name\n"; my $basis2 = "\$basis2\ntype=general\ndefault=$basis_name\n\nset=$basis_name\n"; ELEMENT: foreach my $element_number (eval $atom) { my $charge = $element_number % 2; $charge = 1-$charge if $gen eq 'orca'; my $element = $ATOM[$element_number]; print "$element "; if ($gen eq 'g09') { open L, '>', "$element.com" or die "Can't wtite $element.com: $!\n"; print L "# test BLYP/$basis_name DenFit GFInput Guess=Only Pop=None\n\n$element\n\n$charge 1\n$element\n\n\n"; close L; system($g09, "$element.com"); # my @args = ($g09, "$element.com"); # system(@args) == 0 or die "\nSystem\n@args\nfailed: $?\n"; open L, '<', "$element.log" or die "Can't open $element.log: $!\n"; my (@g09_basis, @g09_fit_basis); while (<L>) { if (m/ Atomic number out of range/) { print "\b"x(length($element)+1); unlink "$element.com", "$element.log" unless $no_del; next ELEMENT; } if (m/ No stored fitting set matched to this AO basis/) { print "No stored fitting set matched to this AO basis\n"; unlink "$element.com", "$element.log" unless $no_del; last ELEMENT; } if (m/ AO basis set in the form of general basis input/ .. m/ \*\*\*\*/) { push @g09_basis, $_; } if (m/ Standard density basis:\s+(\S+)/) { $info .= "\t\t\t\tbasis2: $1\n"; } if (m/ Density basis set in the form of general basis input/ .. m/ \*\*\*\*/) { push @g09_fit_basis, $_; } } close L; if (@g09_basis > 2) { shift @g09_basis; pop @g09_basis; } else { print STDERR "AO basis set $basis_name for $element is wrong\n"; next; } if (@g09_fit_basis > 2) { shift @g09_fit_basis; pop @g09_fit_basis; } else { print STDERR "Density basis $basis_name for $element set is wrong\n"; next; } unlink "$element.com", "$element.log" unless $no_del; $info .= "$element\t "; my @b = g2pribas(@g09_basis); my $sum = 0; foreach (@b) { $sum += @$_; } $basis .= sprintf "%9d%4d\n", $element_number, $sum; for (my $i=0; $i<@b; $i++) { foreach (@{$b[$i]}) { $basis .= sprintf "%d%3d\n", $i, 0+@$_; $info .= 0+@$_; foreach (@$_) { $basis .= join(' ', map {Lprintf($_)} @$_) . "\n"; } } $info .= '/'; } $info =~ s/\/$/\n/; #pp @b; my @b2 = g2pribas(@g09_fit_basis); my $sum2 = 0; foreach (@b2) { $sum2 += @$_; } $basis2 .= sprintf "%9d%4d\n", $element_number, $sum2; for (my $i=0; $i<@b2; $i++) { foreach (@{$b2[$i]}) { $basis2 .= sprintf "%d%3d\n", $i, 0+@$_; foreach (@$_) { $basis2 .= join(' ', map {Lprintf($_)} @$_) . "\n"; } } } # print "@g09_fit_basis\n"; #pp @b2; } elsif ($gen eq 'orca') { open L, '>', "$element.inp" or die "Can't wtite $element.inp: $!\n"; print L "! PBE $basis_name PrintBasis NoIter MiniPrint\n\n* xyz $charge 1\n $element 0 0 0\n*\n"; close L; open F, '-|', "$orca $element.inp 2>/dev/null" or die "Can't run $orca: $!\n"; my (@b,@b2,$inf,$inf2); my $count = 0; my $count2 = 0; while (<F>) { last if m/One Electron integrals/; if (m/There is no basis function on atom/) { print "\b"x(length($element)+1); unlink glob("$element.*") unless $no_del; close F; next ELEMENT; } if (m/(\S+) contracted to (\S+) pattern/) { if (! $inf) { $inf = "{$1}/[$2]"; } else { $inf2 = "{$1}/[$2]"; } } if (m/^ # Basis set for element/..m/end;/) { if (m/^ ([D-S]) (\d+)/) { push @b, sprintf "%d%3d\n", $SPD{$1}, $2; $count++; for (my $i=1; $i<=$2; $i++) { my ($alpha,$d) = map {Lprintf($_)} (split ' ', <F>)[1,2]; push @b, "$alpha $d\n"; } } } if (m/^ # Auxiliary basis set for element/..m/end;/) { if (m/^ ([D-S]) (\d+)/) { push @b2, sprintf "%d%3d\n", $SPD{$1}, $2; $count2++; for (my $i=1; $i<=$2; $i++) { my ($alpha,$d) = map {Lprintf($_)} (split ' ', <F>)[1,2]; push @b2, "$alpha $d\n"; } } } } close F; unless (@b) { print STDERR "No basis\n"; unlink glob("$element.*") unless $no_del; next BASIS if $element eq 'H'; next ELEMENT; } unless (@b2) { print STDERR "No auxiliary basis\n"; next ELEMENT; } unshift @b, sprintf "%9d%4d\n", $element_number, $count; unshift @b2, sprintf "%9d%4d\n", $element_number, $count2; $basis .= join '', @b; $basis2 .= join '', @b2; $info .= sprintf(" %-5s%-30s%s\n",$element,$inf,$inf2) if $inf; unlink glob("$element.*") unless $no_del; } } my %h; $info = join "\n", grep {! $h{$_}++} split "\n", $info if $gen eq 'g09'; open L, '>', $basis_name or die "Can't open $basis_name: $!\n"; print L "$info\n\n"; print L "$basis\n\$end\n"; print L "$basis2\n\$end\n"; close L; print "\n"; } sub g2pribas { my @pribas; while (@_) { $_ = shift @_; if (m/^\s*[A-Z]/) { my ($shells, $contractations, $scale_factor, $unknown) = split; if ($scale_factor != 1 or $unknown != 0) { print; return undef; } my (@alpha,@d); for (my $i=0; $i<$contractations; $i++) { ($alpha[$i],@{$d[$i]}) = split ' ', shift @_; #print "$alpha[$i],@{$d[$i]}\n"; if (length($shells) != @{$d[$i]}) { print STDERR "$shells but @{$d[$i]}\n"; return undef; } } #pp @d; my $count=0; foreach (split '', $shells) { my @bas; for (my $i=0; $i<$contractations; $i++) { $bas[$i] = [$alpha[$i], $d[$i][$count]]; } #pp @bas; push @{$pribas[$SPD{$_}]}, \@bas; $count++; } } } return @pribas; } sub Lprintf { my $num = shift; $num =~ s/d/e/i; my ($m,$e) = split 'e', sprintf("%+.12e", $num); $m =~ s/\.//; $m =~ s/^([+-])/$1./; $e = sprintf "%+03d", $e+1; return $m.'00000e'.$e; }