Новосибирский институт органической химии им. Н.Н. Ворожцова СО РАН
Лаборатория изучения механизмов органических реакций
 |
 |
 |
 |
 |
 |
pri2mol
#!/usr/bin/perl -ws our($asymm,$debug,$all); $debug ||= 0; $all ||= 0; #use Data::Dump 'pp'; # Если нет параметров, то печатаем справку if (! @ARGV or $ARGV[0] eq '-h') { (my $program = $0) =~ s/.*[\/\\]//; print "\n Usage: $program file [files]\n Зависимости: perl\n Этот перл-скрипт читает выдачи квантово-химической программы ПРИРОДА (автор Д.Н. Лайков) и делает файлы, понятные molden'у и заготовки для следующего расчета. Проверенные версии ПРИРОДЫ: p210, p407, p6, p11 (linux), p202 (Windows).\n Создаются следующие файлы: file.xyz декартовы координаты с зарядами для всех шагов оптимизации или для каждой точки irc; file.mos собственные вектора (орбитали) в формате molden'а (почему-то обрабатываются molden'ом очень медленно); file.freq колебания в формате molden'а (при task=hessian); file.MOL оптимизированная геометрия и гессиан, которые можно вставить в след. расчет (не для molden'а); file.ppm хим. сдвиги (sigma и, если посчитан стандарт, то и delta) (не для molden'а); file.xyzppm геометрия для molden'а, с хим. сдвигами (delta) вместо зарядов (molden желательно немного переделать, чтобы он понимал трехзначные \"заряды\" (знаю, как)).\n\n"; exit; } # Вгоняем таблицу Менделеева @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 ); %HoF = ( PBE => {L1=>{C=>-38.0843685585, H=>-0.5872542349, N=>-54.7271355293, O=>-75.1115133857}, }, ); my $num = qr/-?\d+\.\d+/; NEW_FILE: foreach $file (@ARGV) { my (@xXX,@yXX,@zXX); my %formula; my (@thermo,@intensity); my @crit_points; undef $hessian; # Открываем файлы $ppm = "$file.ppm"; $mos = "$file.mos"; $freq = "$file.freq"; $mol = "$file.MOL"; $xyz = "$file.xyz"; $xyzppm = "$file.xyzppm"; open (INP, $file) || die "Cannot open input file $file"; open (MOS, ">$mos") || die "Cannot open input file $mos"; open (FREQ, ">$freq") || die "Cannot open input file $freq"; open (MOL, ">$mol") || die "Cannot open input file $mol"; open (XYZ, ">$xyz") || die "Cannot open input file $xyz"; open (PPM, ">$ppm") || die "Cannot open input file $ppm"; open (XYZPPM, ">$xyzppm") || die "Cannot open input file $xyzppm"; # по дефолту вид расчета и # ограничитель, после нахождения которого нужно печатать геометрию # (для IRC и SCAN печатаем только оптимизированные на каждом шаге геометрии) $task = 'optimize'; $delimiter = 'eng>\$end'; # scan input file my $opt_conv = 0; while (<INP>) { # program Priroda version 4.07 (14.09.2004) # program Priroda version 2.10 (14.01.2002) # Priroda 6 (2006.08.20) if (m/\bPriroda\b\D*(\d+)/) { $pri_ver = $1; } if (/THERMOSTAT/) { close MOS; close FREQ; close MOL; close PPM; close XYZPPM; unlink($mos,$freq,$mol,$ppm,$xyzppm); seek (INP, 0, 0); thermostat(); close INP; close XYZ; next NEW_FILE; } if (/Intrinsic Reaction Coordinate/) {$task = 'irc'}; if (/SCAN\>/ | /task=scan/) {$task = 'scan'}; if (/OPTIMIZATION CONVERGED/) {$opt_conv++}; if (/NMR Shielding Tensors/) {$task = 'nmr'}; if (/COMPUTATION OF SECOND DERIVATIVES/) {$task = 'hessian'}; #print "$task\n"; #if (/Calculation of NMR Chemical Shifts/) {$task = 'nmr'}; # Считываем название основного базиса $basis = $1 if m/^ Basis set input: '(\S+)'/i; $basis =~ s/.*\/// if $basis; # Дополнительный базис $basis2 = $1 if m/^ Basis set input: '(\S+)'/i; $basis2 =~ s/.*\/// if $basis; # Считываем число атомов $natoms = $1 if m/^\s*(\d+)\s+atoms,\s*\d+\s+electrons/; $natoms = $1 if m/Number of atoms\s+(\d+)/; # Заряд молекулы и ее мультиплетность $Charge = $1 if /Charge of molecule\s+(-?\d+)/; $Mult = $1 if /Spin multiplicity\s+(\d+)/; # if (/molecule input:/) { # $_=<INP>; # ($natoms, undef)=split; # } # Функционал $functional = $1 if m/Approximation to E\(xc\)\s+(\S+)/; if (m/formula: (\w+)/) { my $brutto = $1; $brutto =~ s/([A-Z][a-z]?)(?!\d)/${1}1/g; %formula = split /(\d+)/, $brutto; } $ZPVE = $1 if /ZPVE =\s*(\S+)/; if (m/^\s*T =\s*(\S+)\s+K/) { my $T = $1; $T =~ s/0+$//; $T =~ s/\.$//; push @thermo, {T => $T}; <INP>;<INP>;<INP>;<INP>;<INP>; $_ = <INP>; my @a = split; $_ = sprintf("%.2f",$_) for @a[3..5]; $thermo[-1]{lnQ} = $a[1]; $thermo[-1]{S} = $a[2]; $thermo[-1]{U} = $a[3]; $thermo[-1]{H} = $a[4]; $thermo[-1]{G} = $a[5]; } if (m/qm file:/) { $qm = 1; $delimiter = 'mol>\$end'; } } if ($task ne 'irc' && $opt_conv > 1) { $task = 'scan' } if ($task eq 'irc' || $task eq 'scan') { $delimiter = 'OPTIMIZATION CONVERGED'; } if ($task eq 'nmr') { $delimiter = 'Atomic Coordinates:'; } if ($all and $task eq 'irc' || $task eq 'scan') { $delimiter = 'eng>\$end'; $delimiter = 'mol>\$end' if $qm; } print "\$natoms = $natoms\n" if ($debug >=1); print "\$task = $task\n" if ($debug >=1); print "\$theory = QM\n" if ($debug >=1 and $qm); print "\$delimiter = $delimiter\n" if ($debug >=1); # Переходим в начало файла seek (INP, 0, 0); # (флаг, найдены ли заряды и спины) $mulOK=0; # Вновь сканируем файл while (<INP>) { # Читаем координаты if (/Atomic Coordinates:/) { for ($i=1; $i<=$natoms; $i++) { $_=<INP>; ($charge[$i], $x[$i], $y[$i], $z[$i]) = split; } } #print "@charge\n" if ($debug >=1); # Считываем энергию $energy=$1 if $qm && /^ Energy:\s*(\S+)/; $energy=0+$1 if m/^ E =\s*(\S+)/; $energy=0+$1 if m/^eng> E=(\S+)/; $dipole = $1 if m/ dipole = .+? (\S+) D/; if ($task eq 'irc') { if (m/^IRC> point.+? g = (\d+\.\d+)/) { $gradient = $1; } } # Читаем заряды и спины if (m/atom charge spin\s*/) { $mulOK=1; for ($i=1; $i<=$natoms; $i++) { $_=<INP>; (undef, undef, $mulcharge[$i], $spin[$i]) = split; } } if (m/Critical Points of the Density/) { my @tmp; while (<INP>) { last if m/Atomic Coordinates:/; if (m/r=\s*($num),\s*($num),\s*($num)\s*f=($num)/) { push @tmp, ['XX', $1,$2,$3,$4]; } } push @crit_points, \@tmp; } if (m/Geometry Optimization Step (\d+)/) { if (defined($step) && $step > $1) { #print "$step > $1\n"; $new_optimization = 1; } $step = $1; } $geo_opt = 1 if m/^ OPTIMIZATION CONVERGED/; # Печатаем число атомов, энергию и координаты # (а также зарядамы и спины, если они найдены) # if (/$delimiter/ or ($new_optimization and !$geo_opt)) { if (/$delimiter/) { #print "($new_optimization and !$geo_opt)\n"; $new_optimization = 0; $geo_opt = 0; printf XYZ " %d\n", @crit_points ? $natoms+@{$crit_points[-1]} : $natoms; print XYZ " Energy $energy"; print XYZ " Charge $Charge" if $Charge; print XYZ " Mult $Mult" if $Mult && $Mult != 1; print XYZ " Dipole $dipole" if $dipole; undef $dipole; if ($task eq 'irc' and ! $all and $gradient) { print XYZ " Gradient $gradient"; } if ($ZPVE) { print XYZ " ZPE $ZPVE"; } if (@thermo) { foreach (@thermo) { print XYZ " G($_->{T}) $_->{G}"; } } if ($energy && $functional) { if (! %formula) { foreach my $at (@charge) { next unless $at; #print "$at\n"; $formula{$ATOM[$at]}++; } } my $HoF = $energy; foreach my $element (keys %formula) { if (! exists $HoF{$functional}{$basis}{$element}) { undef $HoF; last; } $HoF -= $HoF{$functional}{$basis}{$element}*$formula{$element}; } if (defined $HoF) { printf XYZ " HoF %.2f kcal", $HoF*627.5; } } print XYZ "\n"; #print XYZ ($mulOK ? " Charge Spin\n" : "\n"); if ($mulOK) { for ($i=1; $i<=$natoms; $i++) { printf XYZ "%2.4s %12.8f %12.8f %12.8f %10.4f %8.4f\n", $ATOM[$charge[$i]], $x[$i], $y[$i], $z[$i], $mulcharge[$i], $spin[$i]; } } else { for ($i=1; $i<=$natoms; $i++) { printf XYZ "%2.4s %12.8f %12.8f %12.8f\n", $ATOM[$charge[$i]], $x[$i], $y[$i], $z[$i]; } } if (@crit_points) { my $aref = pop @crit_points; #pp $aref; for ($i=0; $i<@$aref; $i++) { printf XYZ "%2.4s %12.8f %12.8f %12.8f %10.4f\n", @{$aref->[$i]}; } } } # При простой оптимизации геометрии (в т.ч. saddle) if ($task eq 'optimize') { print XYZ if (/OPTIMIZATION CONVERGED/); } # Создаем файлы орбиталей и частот колебаний и геометрия+гессиан print MOS if s/mos>(.*$)/$1/; print MOL if s/MOL>(.*$)/$1/; $hessian=1 if (/SECOND DERIVATIVES/); if ($hessian) { print MOL if s/eng>(.*$)/$1/; } $vibration=1 if (/^ VIBRATIONAL ANALYSIS and THERMOCHEMISTRY/); if ($vibration && $pri_ver) { if ($pri_ver == 2) { push @intensity, split if s/^ ir i\.//; } # if ($pri_ver == 4 or $pri_ver == 6) { if ($pri_ver >= 4) { if (/^ \| Mode \| Freq\. \| Mass\. \| IR Int\./) { <INP>; while (($_=<INP>) =~ /^\s*\|\s*\d+\s*\|/) { push @intensity, (split /\|/,$_)[4]; } } } print FREQ if s/freq>//; } # Разбираемся с хим. сдвигами if (/Calculation of NMR Chemical Shifts/) { <INP>; if (<INP> =~ /additional points/) { while (my($x,$y,$z) = split ' ', <INP>) { push @xXX, $x; push @yXX, $y; push @zXX, $z; } } } if (/NMR Shielding Tensors/) { %TMS = ( 'PBE' => { 'iglo.bas' =>{'C'=>181.7934,'H'=>31.567, 'Si'=>332.2421}, 'II-ccpv3z.bas'=>{'C'=>181.7934,'H'=>31.567, 'Si'=>332.2421}, 'iglo3z.bas' =>{'C'=>182.088, 'H'=>31.599, 'Si'=>334.1568}, 'II-3z.bas' =>{'C'=>182.088, 'H'=>31.599, 'Si'=>334.1568}, 'L2.bas' =>{'C'=>182.94-7.16, 'H'=>31.179-0.031, 'Si'=>352.9499}, 'L1.bas' =>{'C'=>182.94-7.16, 'H'=>31.179-0.031, 'Si'=>352.9499}, 'L22.bas' =>{'C'=>182.94-10, 'H'=>31.179-0.031, # TMS-empir. 'Si'=>352.9499, # TMS # 'N'=>263.2636, # NH3 L1 geom. 'N'=>228.74, # THEOCHEM 'O'=>186.57+77.3, # PhOH [OMR 1993, 31, 489-494] }, '3z.bas' =>{'C'=>182.000, 'H'=>31.392, 'Si'=>328.8347, 'F'=>318.91, 'N'=>-134.594,'O'=>327.9538, 'Al'=>555.98, 'Ti'=>-1085, 'Cl'=>-20.2, #'B'=> 94.05, # sigma(BF3-Et2O) 'B'=>[-0.935433,92.888], 'P'=>277.08, # sigma(H3PO4) 'P'=>353.96+(-62), # sigma(PMe3)+delta(PMe3) }, '3z_sbk.bas' =>{'C'=>182.000, 'H'=>31.392, 'Si'=>328.8347, 'F'=>318.91, 'N'=>-134.594,'O'=>327.9538, 'Al'=>555.98, 'Ti'=>-1085, 'Cl'=>-20.2, 'B'=> 93.90, # sigma(BF3-Et2O) 'P'=>277.08, # sigma(H3PO4) 'P'=>353.96+(-62), # sigma(PMe3)+delta(PMe3) }, 'ccpv3z.bas' =>{'C'=>183.089, 'H'=>31.199, 'Si'=>335.1516}, 'ccpv4z.bas' =>{'C'=>180.54, 'H'=>31.116, 'Si'=>344.3229}, 'III-3z.bas' =>{'C'=>188.0585+2.1, 'H'=>31.288+0.233}, 'IV-3z.bas' =>{'C'=>188.1792+2.1, 'H'=>31.1894+0.233} }, 'B3LYP' => { '3z.bas' => {'C'=>166.214, 'H'=>29.602, 'Si'=>289.219} } ); foreach my $f (keys %TMS) { foreach my $b (keys %{$TMS{$f}}) { if ($b =~ /(.+)\.bas$/) { $TMS{$f}{$1} = $TMS{$f}{$b}; } } } # print header $basis =~ s/.*[\/\\]//; $basis2 =~ s/.*[\/\\]//; print PPM " basis set '$basis/$basis2' \n"; print PPM " Atom SIGMA DELTA\n"; print XYZPPM ' ',$natoms+@xXX,"\n"; print XYZPPM " Energy $energy\n" if $energy; # Печатаем хим. сдвиги for ($i=1; $i<=$natoms; $i++) { $_=<INP>; # читаем строку #(undef, $number, $atom, undef, undef, $sigma) = split; ($number,$atom,$sigma,$principal) = m/\s*Atom\s*(\d+)\s+(\w+).+?(-?\d+\.\d+)(.+)/; if ($asymm) { # Journal of Molecular Structure: THEOCHEM 953 (2010) 83-90 my ($s11,$s22,$s33) = (split ' ', $principal)[1..3]; # Int. J. Mol. Sci. 2002, 3, ~890 ($s33,$s11,$s22) = sort {abs($b-$sigma)<=>abs($a-$sigma)} ($s11,$s22,$s33); # после сортировки нмчего не изменилось #print "$s33,$s11,$s22 - $sigma\n"; $sigma += ($s22-$s11)/($s33-$sigma); } $atom =~ s/(.+):/$1/; # удаляем двоеточие в названии атома printf PPM "%4d %-4.4s %8.4f", $number, $atom, $sigma; printf XYZPPM "%2.4s %12.8f %12.8f %12.8f", $atom, $x[$i], $y[$i], $z[$i]; if (exists $TMS{$functional}{$basis}{$atom}) { my $tms = $TMS{$functional}{$basis}{$atom}; my $ppm; if (ref($tms) eq 'ARRAY') { $ppm = $sigma*$tms->[0]+$tms->[1]; } else { $ppm = $tms - $sigma; } printf PPM " %10.4f", $ppm; printf XYZPPM "%10.4f", $ppm; } printf PPM "\n"; printf XYZPPM "\n"; <INP>; <INP>; <INP>; #Пропускаем три строки } # Печатаем NICS for ($i=0; $i<@xXX; $i++) { $_=<INP>; # читаем строку ($sigma) = /isotropic=\s+(\S+)/; printf PPM "%4d %-4.4s %8.4f %10.4f\n", $natoms+$i+1, 'XX', $sigma, -$sigma; printf XYZPPM "%2.4s %12.8f %12.8f %12.8f %10.4f\n", 'XX', $xXX[$i], $yXX[$i], $zXX[$i], -$sigma; <INP>; <INP>; <INP>; #Пропускаем три строки } } } # Если не закрыть файлы, то они могут не (сразу) записаться на диск close XYZ; close MOS; close FREQ; close MOL; close PPM; close XYZPPM; # Пустые файлы удаляем if (-z $xyz) {unlink $xyz} if (-z $ppm) {unlink $ppm} if (-z $xyzppm) {unlink $xyzppm} if (-z $freq) {unlink $freq} if (-z $mos) {unlink $mos} if (-z $mol) {unlink $mol} # Дописываем интенсивности ИК if (@intensity) { if (-f $freq) { open (FREQ,">>$freq") || die "Cannot open $freq to add"; print FREQ "[INT]\n"; foreach (@intensity) { print FREQ " $_\n"; } close FREQ; } } # if (@crit_points) { # pp @crit_points; # } } sub thermostat { my (@mols,@trj); while (<INP>) { if (/^ point /) { my @t = /(-?\d+(?:\.\d+)?)/g; push @trj, \@t; } if (/Atomic Coordinates:/) { my @mol; $mol[0] = $trj[-1] if @trj; while (<INP>) { last if /#/; s/\s*(\d+)\s/ $ATOM[$1] /; push @mol, $_; } push @mols, \@mol; } } unshift @{$mols[0]}, $trj[0]; my $N = $#{$mols[0]}; foreach my $mol (@mols) { print XYZ "$N\n"; print XYZ " Energy $mol->[0][2] point $mol->[0][0] s $mol->[0][1] g $mol->[0][3]\n"; shift @$mol; print XYZ @$mol; } open TRJ, '>', "$file.trj" or die "Can't open $file.trj: $!\n"; foreach (@trj) { printf TRJ "%4d %10f %14.8f %12.8f\n", $_->[0], $_->[1], $_->[2], $_->[3]; } close TRJ; }