Новосибирский институт органической химии им. Н.Н. Ворожцова СО РАН
Лаборатория изучения механизмов органических реакций
 |
 |
 |
 |
 |
 |
conformers
#!/usr/bin/perl -ws use strict; #use Data::Dump; #use List::Util 'sum'; ################################################################################## # Parameters for mopac # MOPAC my $mopac_x = '/opt/mopac/MOPAC2016.exe'; $ENV{'MOPAC_LICENSE'} = '/opt/mopac/'; my $mop_key = 'XYZ PRECISE EF'; my $cxcalc_x = '/usr/local/lib/MarvinBeans/bin/cxcalc'; $ENV{'CHEMAXON_LICENSE_URL'} = "$ENV{'HOME'}/.chemaxon/license.cxl"; my $maxconformers = 80; my $diversity = 0.2; my $timelimit = 36000; my $optimization = 2; #my $optimization = 1; # normal (limit=0.001) ## Parameters for Priroda 11 #my $p11_x = '/usr/local/lib/priroda/p1'; #my $p11_basis = '/usr/local/lib/priroda/basis/qm.in'; # Parameters for Priroda 14 my $mpiexec_x = '/usr/local/lib/priroda/mpiexec'; my $p14_x = '/usr/local/lib/priroda/p'; my $p14_basis = '/usr/local/lib/priroda/basis/qm.in'; # Parameters for tinker scan. my $tinker_dir = '/usr/local/lib/tinker'; #my $scan_x = "$tinker_dir/source/scan.x"; my $scan_x = "$tinker_dir/bin/scan"; my $tinker_params = "$tinker_dir/params"; #my $tinker_options = "$tinker_params/mmff 0 5 20 0.0001"; # MMFF94s and default #my $tinker_options = "$tinker_params/mm2 0 5 20 0.0001"; # MM2 and default # Default options for tinker scan our $tinker_rms ||= 0.0001; our $tinker_mm ||= 'mmff'; our $tinker_tors ||= 0; our $tinker_directions ||= 5; our $tinker_energy ||= 20; our $tinker_rotate ||= 0; # параметры (числа через пробел): # 0 - Selection of Torsional Angles for Rotation (0 - Automatic) # 5 - the Number Search Directions for Local Search # 20 - the Energy Threshold for Local # 0.0001 - RMS Gradient per Atom Criterion # $tinker_mm eq 'mm2' : babel needs # $tinker_mm eq 'mmff' : sdf2tinkerxyz needs #http://sdf2xyz2sdf.sourceforge.net/?Description #my $sdf2xyz2sdf_x = "$ENV{'HOME'}/Programs/open3dtools/bin/sdf2tinkerxyz"; my $sdf2xyz2sdf_x = "$tinker_dir/bin/sdf2tinkerxyz"; my $tminimize_x = "$tinker_dir/bin/newton"; #my $tminimize_x = "$tinker_dir/bin/optimize"; #my $tminimize_x = "$tinker_dir/bin/minimize"; # Parameters for vconf. my $vconf_x = "$ENV{'HOME'}/Programs/vconf_v2_acad/vconf.exe"; my $vconf_options = ''; #'-sr b' Stereochemistry Restrictions # b Filtering based upon both chirality and cis/trans # n No chirality or cis/trans filtering. # c Filtering based upon only chirality # i Filter based upon only double bond cis/trans # #'-ns 100' numSteps, the chief determinant of the thoroughness # and duration of a calculation in search mode # #'-sw 0.5' searchWidth, between 0 and 1.0 # Больше - шире охватывать возможные конформеры # Меньше - тщательнее исследовать самые стабильные # Parameters for confab. my $confab_x = "/usr/local/bin/obabel"; my $confab_options = ''; # -r <rmsd> RMSD cutoff (default 0.5 Angstrom) # -e <energy> Energy cutoff (default 50.0 kcal/mol) # -c <#confs> Max number of conformers to test (default is 1 million) # -a Include the input conformation as the first conformer # Parameters for babel. my $babel_x = '/usr/local/bin/babel'; my $obminimize_x = '/usr/local/bin/obminimize'; our $obminimize_n ||= 10000; # Parameters for symmetry. my $symmetry_x = '/usr/local/bin/symmetry'; ################################################################################## # Переменные опций our ($logfile,$h,$help,$charge,$mult,$mopac_key,$keep_out,$mmff_type,$method, $tinker,$cxcalc,$confab,$no_confab_energy,$energy,$no_mopac,$vconf,$RMS,$sort_RMS, $no_ret,$only_gen,$no_cf,$MM_key,$ellips, $np, $no_sort_energy,$AW,$QM_SP, ); $method ||= 'RM1'; # Вроде, лишние переменные: $predefined_method (my $program = $0) =~ s/^.*[\/\\]//; if ($h || $help) { print <<H; Анализ конформеров. Данная программа берет молекулярно-механический набор конформеров (или пытается генерироват?? его из единственной структуры, если заданы соответствующие опции), удаляет неуникальные, оставшиеся оптимизирует и снова удаляет неуникальные. Usage: $program [-logfile] [-option=string] file.ext [file1.ext ...] Зависимости: Perl, [MOPAC2009, babel, cxcalc, tinker, vconf]. По умолчанию -- инпут с конформерами в формате xyz, при этом babel и генераторы конформеров (cxcalc, tinker, vconf) не нужны. Если расширение файла не .xyz, то будет попытка конвертации его в xyz с помощью babel. В начале скрипта прописаны пути к используемым внешним программам и параметры вызова последних. Их нужно просмотреть и поменять, если необходимо. Опции: -logfile вывод информации о работе скрипта будет не на консоль, в файл с названием, основанным на имени input-файла и расширением ".log" -no_mopac Вычисления MOPAC'ом не производятся. Только удаление дублей из предоставленного в качестве инпута или из вычисленного молекулярной механикой набора конформеров. Эта опция необходима, например, для объединения уже посчитанных MOPAC'ом результатов после разных генераторов конформеров. -method=RM1|PM6|PM7|AM1|PM3|MNDO Минимизация геометрии будет выполнена mopac'ом По умолчанию -method=RM1 -method=QM Вмест?? MOPAC будет Природа со своим новым полуэмпирическим методом (медленнее, чем MOPAC) -QM_SP Добавлять во 2-ю строку Energy_SP - энергию неоптимизированной структуры -method=MMFF94|UFF|MMFF94s|Ghemical|Gaff Минимизация геометрии будет выполнена молекулярной механикой программой obminimize (openbabel). -obminimize_n default 10000 -method=tinker Минимизация геометрии будет выполнена молекулярной механикой tinker'ом (по умолчанию программа $tminimize_x, силовое поле $tinker_mm, rms $tinker_rms). Эти параметры заданы в начале скрипта. -tinker_rms default 0.00001 -mopac_key='string' Добавочные ключевые слова для mopac. См. help по MOPAC http://openmopac.net/manual/index.html Несколько слов -- через пробел в общих для всей строки кавычках. -keep_out Don't delete inputs and outputs of mopac and Priroda -charge=number Заряд (для mopac и для Природы) По умолчанию charge=0 или 1 (вычисляется с учетом значения -mult) -mult=number Мультиплетность (для mopac и для Природы) -np=number Число ядер (тредов) для mopac и для Природы (-np=1 default) -MM_key='string' Параметры для obminimize. -ff MMFF94, остальные умолчаемые. -energy=number минимальная разница энергий конформеров, которые будут считаться неидентичными. Если эта опция не задана, конформеры будут срав??иваться каждый с каждым. Если энергия в хартри, рекомендуется -energy=0.001 (~0.6 kcal/mol) -ellips=number Среднеквадратичная разница длин осей эллипсоида инерции. Если эта разница у двух сравниваемых структур больше, то структуры считаются считаются различными. Рекомендуется -ellips примерно как -RMS Если эта опция не задана, конформеры будут сравниваться каждый с каждым. -RMS=number минимальное значение RMS, при котором конформеры будут считаться идентичными. RMS = sqrt(sum(m*r**2)/sum(m)). По умолчанию 0.1 -sort_RMS значение RMS, при котором будет проводится перенумерация атомов текущего конформера по предыдущему. Нужно задавать в явном виде, иначе перенумерации не будет. -no_ret не будет сохраняться конфигурация молекулы. Нужно, если имеем дело со смесью оптических изомеров, и их не следует различать. -cxcalc конформеры генерирует консольная утилита MarvinSketch. Нужна академическая лицензия на MarvinSketch (?). См. $cxcalc_x conformers -h Если формат файла такой, в котором есть информация о связях и который понимает MarvinSketch (mol, mol2, mrv, pdb, sdf, syb) то babel не ну??ен. -vconf[='vconf options'] конформеры генерирует Vconf. -vconf='-ns 200' ограничиться ~200 конформерами. Нужна академическая лицензия. Если формат mol, то babel не нужен. -tinker[='parameters'] конформеры генерирует tinker. parameters - это строка вида 'mm2 0 5 100 0.0001' (либо без последних членов, вплоть до 'mm2'), где mm2 - это имя файла в директории $tinker_params. По умолчанию силовое поле mmff (MMFF94), но для этого нужна программа sdf2tinkerxyz из пакета sdf2xyz2sdf и tinker v.6. Если формат txyz (родной формат tinker'а), то babel не нужен. Если tinker надолго задумается, то можно убить скрипт (^C). При этом найденные конформеры останутся в виде файлов. При повторном запуске скрипта (с теми же опциями) tinker уже не будет запускаться, скрипт обработает файлы конформеров и пойдет дальше. -tinker_rot_bonds только зачитать кол-во вращабельных связей и закончить -tinker_tors По умолчанию 0 (автоматически определять вращабельные связи) Если 1, то программа остановится и будет ждать ввода пар атомов (номера через пробел, по паре на строку), специфицирующих вращабельные связи. Конец ввода - пустая строка. Если 2, то то же самое, но для связей, вращение вок??уг которых нужно заморозить (например, С=С). -confab[='confab options'] конформеры генерирует confab. confab не выдает энергий получившихся конформеров (ротамеров). Потому, если не задана опция -no_confab_energy, получившиеся конформеры оптимизируются obminimize (MMFF94). -no_confab_energy confab отрабатывает очень быстро, энергии ротамеров считаются гораздо дольше. Эта опция предотвращает расчет энергий. -mmff_type=a1,t1,a2,t2,... Пары номер_атома,тип_атома При автоматической генерации txyz некоторые mmff-типы могут быть неподходящими. Их можно переопределить. Типы атомов см. в файле $tinker_params/$tinker_mm.prm Например, 165 "PYRIDINIUM-TYPE N"; 128 "H PYRIDINE N(+)" -only_gen Не проводить сравнение конформеров и не запускать mopac. Эта опция уместна для tinker'а, когда конформеров очень много, чтобы получить их в xyz-формате и отсортированные по энергиям. $program -only_gen multi-mol.xyz : просто сортировка. -no_cf Не проводить сравнение конформеров, но считать их и сортировать. -AW Задает другие атомные веса. Нужно, чтобы более точно различать конформеры, отличающиеся только расположением легких атомов. -AW is -AW=H,6 -help То же самое, что и -h плюс некоторые к??мментарии. H } if ($help) { print <<HELP; Эта программа -- не генератор конформеров. Ее предназначение -- пропустить через MOPAC генерированный молекулярной механикой набор конформеров и удалить дубли. Смысл этой процедуры в следующем. Некоторые молекулярно-механические программы (их список см. ниже) умеют проводить конформационный анализ автоматически. Но они генерируют большое кол-во конформеров, многие из которых оказываются одинаковыми (особенно у Marvin), После оптимизации квантовой химией число уникальных конформеров еще сокращается, иногда сильно, иногда не очень. После полуэмпирики довольно высока вероятность того, что в начале будут действительно самые стабильные конформеры, их можно уже считать считать неэмпирикой. Список бесплатных генераторов конформеров. Marvin http://www.chemaxon.com/product/conformer.html Vconf http://www.verachem.com/products/vconf/ TINKER http://dasher.wustl.edu/tinker/ Frog http://bioserv.rpbs.jussieu.fr/cgi-bin/Frog2 Balloon http://users.abo.fi/mivainio/balloon/ OMEGA http://www.eyesopen.com/products/applications/omega.html Multiconf-DOCK http://dock.compbio.ucsf.edu/Contributed_Code/multiconfdock.htm confab http://code.google.com/p/confab/ RDKit http://www.rdkit.org/ Некоторые особенности генераторов ко??формеров, основанные на моем ??пыте: Marvin (плагин Marvin Sketch или cxcalc conformers) имеет сложности с генерированием конформеров соединений циклогексанового ряда с фиксированным положением (R,S) заместителей. Разработчики обещают это исправить; пока рекомендую слежующий алгоритм: найти конформеры нарисованного соединения, сделать из него оптический изомер, найти его конформеры, сделать им reflect, объединить оба набора конформеров. Marvin способен выдавать "конформеры", которые нельзя получить без разрыва связи (например, цис- и транс-декалины), но можно и зафиксировать конфигурацию и получить что-то одно. Используется либо примитивное силовое поле Дрейдинга (нет никаких ограничений на элементы и очень мало - на типы связей), либо MMFF94 (получается вполне реалистично, но не все). У Marvin'а из всех опробованных мной генераторов наиболее широкий охват всевозможных структур (правдоподобность полученных энергий, конечно сомнительна, но она и не очень важна, все равно потом оптимизировать полуэмпирикой). Предварительный набор конформеров проще всего получать именно с помощью этой программы, через графический интерфейс MarvinSketch (Tools - Conformation - Conformers. Это ино??да долго. Нужна академическая лицензия, но она нежесткая, переносится на другие машины.) Tinker (программа scan из его состава) не может генерировать конформеры циклогексана, декалина. У него очень хорошо получается конформационный анализ открыто-цепных и макроциклических соединений. Tinker -- Open Source. Использует честные силовые поля, отсюда и ограничение: как только молекула содержит что-то чуть не стандартное, так ошибка (отсутствие параметров). В этом отношении MMFF94 лучше MM2. Tinker не сохраняет конфигурацию двойной связи и R/S. vconf считает хорошо и очень быстро. Использует силовое поле UFF, "всеядность" его по-меньше, чем у Marvin'а, и конформеров, как правило, выдает по-меньше. Можно зафиксировать конфигурацию при двойной связи и/или у тетраэдрического углерода. Основной недостаток vconf - жесткая академическая лицензия, только на год и привязанная к MAC-адресу. confab - свободный систематический генератор, находит открыто-цепочечные конформеры (ротамеры), но не трогает циклы и гетероатом-Н. Не делает оптимизации геометрии получившихся конформеров (поэтому начальную структуру нужно брать "хорошую"). Конформеров получаетс?? меньше, чем ?? tinker'а. К этим четыем генераторам конформеров в скрипте есть интерфейс. Frog работает через web-интерфейс. Конформеров находит мало. Multiconf-DOCK работает, но конформеров находит мало. OMEGA несмотря на труднопредолимые сложности с получением академической лицензии, конформеров находит мало. Другие генераторы я не пробовал. Программа рассчитана на использование MOPAC2009 http://openmopac.net/Download_MOPAC_Executable_Step2.html Для него нужна академическая лицензия, действующая 1 год, но автоматически продлевающаяся после установки новой версии MOPAC. Главное преимущество MOPAC2009 перед mopac7 - имеются параметризации RM1, в которой мало элементов, но параметры для них очень хорошо подобраны, и PM6, охватывающая чуть ли не всю периодическую таблицу. Эти параметризации являются \"продолжетелями\" AM1 и PM3. Оптимизация геометрии идет в декартовых координатах -- это дольше, чем с z-матрицей, но зато надежнее, что очень важно для автоматической работы. Отсеивание дублей в программе реализовано так: две молекулы ориентируются по осям их тензоров инерции, затем вычисляется RMSD между ними. При этом соответс??вие атомов одной молекулы атомам другой определяется просто по их пространственной близости (заодно делается перенумерация второй молекулы, если она достаточно похожа на первую). Поэтому, если собственные числа двух или всех трех осей тензора инерции одинаковы, то будут ошибки. Такое практически не встречается для несимметричных или низкосимметричных молекул, а вот для высокосимметричных - не исключено. Частично ситуация исправлена, но для некоторых случаев (типа диметилацетилена) проблема остается. Поэтому на всякий случай нужно вручную просматривать получившиеся конформеры и отсеивать лишние. Признак одинаковости - те же самые группы симметрии, очень близкие энергии и эллипсоиды инерции. Еще некоторые замечания. Вместо conformers -logfile file.xyz на linux можно использовать conformers file.xyz | tee file.log Обычно в файле есть энергии конформеров (генераторы их вычисляют). Если конформеров очень много (100 и больше), сравнение их похожести занимает длительное время. Чтобы ускорить процесс, можно задать опцию -energy, тогда конформеры, энергии которых различаются больше чем на заданную величину, сразу будут считаться не??охожими. MOPAC, как и любая ква??товая химия, должен знать правильные заряд и мультиплетность. Если у вас нейтральное соединение или катион в синглетном состоянии, то ничего указывать не нужно, вычислится автоматически из брутто-формулы. В противном случае следует задать заряд опцией -charge и мультиплетность через -mopac_key Параметризация в MOPAC. Лучше всего показала себя RM1. Но в ней очень мало элементов. Поэтому, если параметризация не задана явно через -method, будет RM1, но если есть отсутствующие в ней элементы, то PM6. Результат работы программы -- наборы нумерованных в соответствие с энергией xyz- и arc-файлов. В xyz-файлах во 2-й строке, кроме энергий (в ккал/моль) содержатся еше определенные MOPAC'ом группы симметрии, длины полуосей эллипсоидов инерции, а также номера конформера в первоначальном наборе. Эти номера фигурируют и в output'е программы. Разработка этой программы поддержана РФФИ (грант 13-03-00427a). HELP } $| = 1; $ENV{'GFORTRAN_UNBUFFERED_ALL'} = 1; if ($vconf && $vconf eq '1' && ! $no_ret) { $vconf_options .= '-sr c' unless $vconf_options =~ /-sr/; } my $tinker_MM = $tinker_mm; # Дубль, который не будем модифицировать в дальнейшем my ($calc_charge); # For MOPAC and QM $mult ||= 1; my (%RM1,%PM6,%PM7,%QM, $predefined_method); my @mopac_mult = qw(SINGLET DOUBLET TRIPLET QUARTET QUINTET SEXTET SEPTET OCTET NONET); $mop_key .= " $mopac_mult[$mult-1]" if $mult > 1; $np ||= 1; $mop_key .= " THREADS=$np"; unless ($no_mopac) { if ($method eq 'QM') { if (-x $p14_x) { $predefined_method = 1; $charge ||= 0; %QM = map {$_=>1} qw( H Li Be B C N O F Na Mg Al Si P S Cl Sc Ti V Cr Mn Fe Co Ni Cu ); } else { die "Priroda 14 executable $p14_x not found.\n"; } die "No basis file $p14_basis for Priroda (QM method)\n" if ! -f $p14_basis; } elsif ($method =~ /MMFF94s?|UFF|Ghemical|Gaff/) { if (-x $obminimize_x) { $predefined_method = 1; if ($MM_key) { $MM_key .= " -n $obminimize_n" if $MM_key !~ /-n\b/; } else { $MM_key = "-n $obminimize_n"; } } else { die "obminimize executable $obminimize_x not found.\n"; } } elsif ($method =~ /tinker/) { if (-x $tminimize_x) { $predefined_method = 1; } else { die "tminimize executable $tminimize_x not found.\n"; } } elsif ($method =~ /RM1|PM6|PM7|AM1|PM3|MNDO/) { if (-x $mopac_x) { $predefined_method = 1 if $method; $method ||= 'RM1'; $calc_charge = 1 unless defined $charge; $mopac_key = '' unless defined $mopac_key; if ($mopac_key =~ /\bTS\b/i) { $mop_key =~ s/\bEF\b/TS/i && $mopac_key =~ s/\s*\bTS\b//i; } %RM1 = map {$_=>1} qw(H C N O F P S Cl Br I); %PM6 = map {$_=>1} qw( 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 Lu Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi ); %PM7 = map {$_=>1} qw( 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 Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Ce Pr Nd Sm Eu Gd Tb Dy Ho Er Tm Yb Lu ); } else { print "MOPAC executable $mopac_x not found.\nForce -no_mopac option.\n"; $no_mopac = 1; } } else { print "Unknown method $method. Force -no_mopac option.\n"; $no_mopac = 1; } } if ($energy) { if ($energy =~ /^(\d+(?:\.\d+)?)(kcal|kJ|au)?/i) { if ($2) { $energy = $1*627.5 if lc($2) eq 'au'; $energy = $1/4.186 if lc($2) eq 'kj'; } } else { die "Invalid -energy=$energy\n"; } } $RMS = 0.1 unless defined $RMS; #$sort_RMS = 3 * $RMS unless defined $sort_RMS; # Массы изотопов my %massa = qw( H 1.007947 He 4.0026022 Li 6.9412 Be 9.0121823 B 10.8117 C 12.01078 N 14.00672 O 15.99943 F 18.99840325 Ne 20.17976 Na 22.9897702 Mg 24.30506 Al 26.9815382 Si 28.08553 P 30.9737612 S 32.0655 Cl 35.4532 Ar 39.9481 K 39.09831 Ca 40.0784 Sc 44.9559108 Ti 47.8671 V 50.94151 Cr 51.99616 Mn 54.9380499 Fe 55.8452 Co 58.9332009 Ni 58.69342 Cu 63.5463 Zn 65.4094 Ga 69.7231 Ge 72.641 As 74.921602 Se 78.963 Br 79.9041 Kr 83.7982 Rb 85.46783 Sr 87.621 Y 88.905852 Zr 91.2242 Nb 92.906382 Mo 95.942 Ru 101.072 Rh 102.905502 Pd 106.421 Ag 107.86822 Cd 112.4118 In 114.8183 Sn 118.7107 Sb 121.7601 Te 127.603 I 126.904473 Xe 131.2936 Cs 132.905452 Ba 137.3277 La 138.90552 Ce 140.1161 Pr 140.907652 Nd 144.243 Sm 150.363 Eu 151.9641 Gd 157.253 Tb 158.925342 Dy 162.5001 Ho 164.930322 Er 167.2593 Tm 168.934212 Yb 173.043 Lu 174.9671 Hf 178.492 Ta 180.94791 W 183.841 Re 186.2071 Os 190.233 Ir 192.2173 Pt 195.0782 Au 196.966552 Hg 200.592 Tl 204.38332 Pb 207.21 Bi 208.980382 Th 232.03811 Pa 231.035882 U 238.028913 ); if (defined $AW) { $AW = 'H,6' if $AW eq '1'; my %AW = (split /,| /, $AW); #dd %AW; foreach my $at (keys %AW) { if ($at =~ /\W/ or $AW{$at} !~ /^\d+(\.\d+)?$/ ) { warn "Invalid -AW= argument: $at,$AW{$at}\n"; next; } $massa{$at} = $AW{$at}; } } MOL: foreach (@ARGV) { my $file_name = $_; unless (-f $file_name) { warn "No $file_name\n"; next; } my ($name,$ext) = /(.+)\.(.*)/; my $EXT = $ext; my @mols_mm; my $tinker_options = ''; if ($logfile) { if (open LOG, '>', "$name.log") { select LOG; $| = 1; } else { warn "Can't write to $name.log: $!\nPrint to stdout\n"; } } ########################################################################## if (defined $cxcalc) { die "No cxcalc executable $cxcalc_x\n" unless -x $cxcalc_x; # Наверное, этот кусок не нужен, cxcalc и так понимает все форматы # if (defined $babel && $babel ne 'mol') { # die "No babel executable $babel_x\n" unless -x $babel_x; # my @args = ($babel_x, "-i$babel", "$name.$ext", "-omol", "$name.mol"); # print "\nRun babel to tranform $babel to mol for cxcalc:\n@args\n"; # system(@args) == 0 or do {warn "System\n@args\nfailed: $?\n"; next}; # $ext = 'mol'; # } my @args = ($cxcalc_x, "$name.$ext", '-o', "$name.cxcalc.xyz", '-v', 'conformers', '-f', 'xyz', '-m', $maxconformers, '-d', $diversity, '-l', $timelimit, '-O', $optimization); print "\nRun cxcalc:\n@args\n"; print 'Time limit ', $timelimit/60, " min\n"; system(@args) == 0 or do {warn "System\n@args\nfailed: $?\n"; next}; @mols_mm = read_molden("$name.cxcalc.xyz"); unlink("$name.cxcalc.xyz"); unlink "$name.mol" if $file_name ne "$name.mol"; } ########################################################################## elsif (defined $tinker) { # предыдущие выдачи tinker'а my @tinker_files = glob("$name.[0-9][0-9][0-9] $name.[0-9][0-9][0-9][0-9] $name.[0-9][0-9][0-9][0-9][0-9]"); #print "@tinker_files\n"; if (! @tinker_files) { die "No tinker scan executable $scan_x\n" unless -x $scan_x; #die "No force field $tinker_params for tinker\n" unless -r "$tinker_params.prm"; if ($tinker ne '1') { my ($mm, $tors, $directions, $energy, $rms) = split ' ', $tinker; # Ищем силовое поле тинкера if (-f "$tinker_params/$mm.prm") { $mm = "$tinker_params/$mm"; } elsif (-f $mm) { $mm =~ s/\.prm$// or die "Not prm file $mm\n"; } elsif (-f "$tinker_params/$mm") { $mm =~ s/\.prm$// or die "Not prm file $tinker_params/$mm\n"; $mm = "$tinker_params/$mm"; } elsif (-f "$mm.prm") { $mm = "./$mm"; } else { die "Can't find tinker FF $mm\n"; } die "Can't read FF $mm.prm for TINKER\n" unless -r "$mm.prm"; $tinker_mm = $mm; $tinker_tors = $tors if defined $tors; $tinker_directions = $directions if defined $directions; $tinker_energy = $energy if defined $energy; $tinker_rms = $rms if defined $rms; } else { $tinker_mm = "$tinker_params/$tinker_mm" if $tinker_mm !~ /$tinker_params/; } $tinker_options = "$tinker_mm $tinker_tors $tinker_directions $tinker_energy $tinker_rms"; #die "Invalid parameters for TINKER: $tinker_options\n" if # $tinker_options !~ /^\S+\s+\d+\s+\d+\s+\d+\s+\d+\.\d+/; if ($ext ne 'txyz') { if ($tinker_mm =~ /mm2$/) { warn("No babel executable $babel_x\n"),next unless -x $babel_x; my @args = ($babel_x, '-f', '1', '-l', '1', '-h', "$name.$ext", "$name.txyz"); print "\nRun babel to tranform $ext to txyz for tinker scans:\n@args\n"; system(@args) == 0 or do {warn "System\n@args\nfailed: $?\n"; next}; $ext = 'txyz'; } elsif ($tinker_mm =~ /mmffs?$/) { unless (-x $sdf2xyz2sdf_x) { warn "No sdf2xyz2sdf executable $sdf2xyz2sdf_x to generate txyz fo MMFF94.\n"; warn "Try -tinker=mm2 option.\n"; next; } #warn("No sdf2xyz2sdf executable $sdf2xyz2sdf_x\n"),next unless -x $sdf2xyz2sdf_x; if ($ext ne 'sdf') { warn("No babel executable $babel_x\n"),next unless -x $babel_x; my @args = ($babel_x, '-f', '1', '-l', '1', '-h','--title','', "$name.$ext", "$name.sdf"); print "\nRun babel to tranform $ext to sdf for sdf2tinkerxyz:\n@args\n"; system(@args) == 0 or do {warn "System\n@args\nfailed: $?\n"; next}; } else { # delete comments open L, '<', "$name.sdf" or die "Can't open $name.sdf: $!\n"; my @sdf = <L>; close L; $sdf[0] = $sdf[1] = "\n"; open L, '>', "$name.sdf" or die "Can't write $name.sdf: $!\n"; print L @sdf; close L; } my $arg = "$sdf2xyz2sdf_x < $name.sdf"; print "\nRun sdf2tinkerxyz to tranform sdf to txyz for tinker scans:\n$arg\n"; system($arg) == 0 or do {warn "System\n$arg\nfailed: $?\n"; next}; rename 'NONAME.xyz', "$name.txyz"; rename 'NONAME.key', "$name.key"; if (! $no_ret) { open L, '>>', "$name.key" or die "Can't open $name.key: $!\n"; print L "ENFORCE-CHIRALITY\n"; close L; } } else { warn "$name.$ext may be transformed to txyz only for mm2 and mmff FF.\n"; next; } } if (defined $mmff_type) { my %atom_type = split ',', $mmff_type; open TXYZ, '<', "$name.txyz" or do {warn "Can't open $name.txyz: $!\n"; next}; open TMP, '>', "$name.t___" or do {warn "Can't write to file: $!\n"; next}; while (<TXYZ>) { while (my ($atom,$type) = each %atom_type) { if (s/^(\s*$atom\s+\S+\s+\S+\s+\S+\s+\S+\s+)\d+/${1}$type/) { delete $atom_type{$atom}; last; } } #print; print TMP; } close TMP; close TXYZ; rename "$name.t___", "$name.txyz"; } my $rot_bonds_file = 'rot_bonds_file'; if ($tinker_rotate) { $tinker_options =~ s/^\s*(\S+).*/$1/; #$tinker_rotate =~ s/\D+/ /g; open L, '>', $rot_bonds_file or die "Can't open $rot_bonds_file: $!\n"; print L "1\n"; while ($tinker_rotate =~ m/(\d+)\D+(\d+)/g) { print L $1<$2 ? "$1 $2" : "$2 $1", "\n"; } print L "\n$tinker_directions\n$tinker_energy\n$tinker_rms\n"; close L; } my @args = ($scan_x, "$name.txyz", split(' ', $tinker_options)); print "\nRun tinker scan:\n@args\n"; if (-f $rot_bonds_file) { open TINKER, "@args < $rot_bonds_file |" or die "Can't run @args < $rot_bonds_file : $!\n"; } else { open TINKER, '-|', "@args" or do {warn "Can't run @args : $!\n"; next}; } $| = 1; #select(((select(TINKER), $|=1))[0]); my @minima; while (<TINKER>) { if (/Torsions/) { print; } print if /^ INITROT/; if (/Potential Surface Map\s+Minimum\s+(\d+)\s+(\S+)/) { print; open TE, '>>', "$name.TE" or die "Can't open $name.TE: $!\n"; print TE "$1 \t$2\n"; close TE; push @minima, [$1,$2]; } } close TINKER; for (my $count=0; $count<@minima; $count++) { #my $file = $name . '.' . sprintf("%04d", $count+1); my $file = $name . '.' . sprintf("%03d", $minima[$count][0]); my $Energy = $minima[$count][1]; my $mol = txyz2xyz($file,$Energy); push @mols_mm, $mol; } } else { open TE, '<', "$name.TE" or die "Can't open $name.TE: $!\n"; foreach my $file (@tinker_files) { my $mol; my $l = <TE>; chomp $l; ($mol->[0]{'Energy'}) = (split ' ', $l)[1]; open F, '<', $file or do {warn "Can't open $file: $!\n"; next}; <F>; while (<F>) { my ($at,$x,$y,$z) = (split)[1..4]; chop $at until element($at); push @$mol, [$at,$x,$y,$z]; } close F; push @mols_mm, $mol; } close TE; } unlink "$name.txyz" if $file_name ne "$name.txyz"; unlink "$name.sdf" if $file_name ne "$name.sdf"; unlink glob("$name.[0-9][0-9][0-9] $name.[0-9][0-9][0-9][0-9] $name.[0-9][0-9][0-9][0-9][0-9]"); unlink "$name.TE", "$name.key"; #write_molden(@mols_mm, "${name}_tinker.xyz"); } ########################################################################## elsif (defined $vconf) { die "No vconf executable $vconf_x\n" unless -x $vconf_x; if ($ext ne 'mol') { warn("No babel executable $babel_x\n"),next unless -x $babel_x; my @args = ($babel_x, '-f', '1', '-l', '1', "$name.$ext", "$name.mol", '--title', ''); print "\nRun babel to tranform $ext to mol for Vconf:\n@args\n"; system(@args) == 0 or do {warn "System\n@args\nfailed: $?\n"; next}; $ext = 'mol'; } $vconf = '' if $vconf eq '1'; $vconf_options .= " $vconf"; my @args = ($vconf_x, split(' ', $vconf_options), "$name.$ext"); print "\nRun vconf:\n@args\n"; system(@args) == 0 or do {warn "System\n@args\nfailed: $?\n"; next}; (my $vconf_name = $name) =~ s/\..*//; #my $sdf = -f "${name}_mol_1_confs.sdf" ? "${name}_mol_1_confs.sdf" : "${name}_vconf.sdf"; open SDF, '<', $vconf_name.'_mol_1_confs.sdf' or die "Can't open ${vconf_name}_mol_1_confs.sdf: $!\n"; my $count = 0; while (<SDF>) { if (m/^\s*(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(\w+)/) { push @{$mols_mm[$count]}, [$4,$1,$2,$3]; next; } if (m/^> <Energy>$/) { if (<SDF> =~ /^^\s*(-?\d+\.\d+)$/) { unshift @{$mols_mm[$count]}, {Energy => $1}; next; } } if (m/^\$\$\$\$$/) { unshift @{$mols_mm[$count]}, {} if ! exists $mols_mm[$count][0]{'Energy'}; $count++; } } close SDF; unlink("$name.log", "${vconf_name}_vconf.sdf", "${vconf_name}_mol_1_confs.sdf"); unlink "$name.mol" if $file_name ne "$name.mol"; #write_molden(@mols_mm, "${name}_vconf.xyz"); } ########################################################################## elsif (defined $confab) { die "No confab executable $confab_x\n" unless -x $confab_x; if ($ext ne 'mol') { warn("No babel executable $babel_x\n"),next unless -x $babel_x; my @args = ($babel_x, '-f', '1', '-l', '1', "$name.$ext", "$name.mol", '--title', ''); print "\nRun babel to tranform $ext to mol for confab:\n@args\n"; system(@args) == 0 or do {warn "System\n@args\nfailed: $?\n"; next}; $ext = 'mol'; } $confab = '' if $confab eq '1'; $confab_options .= " $confab"; my @args = ($confab_x, "$name.$ext", '-O', "${name}_confab.sdf", '--confab', split(' ', $confab_options)); print "\nRun confab:\n@args\n"; system(@args) == 0 or do {warn "System\n@args\nfailed: $?\n"; next}; open SDF, '<', "${name}_confab.sdf" or die "Can't open ${name}_confab.sdf: $!\n"; my $count = 0; while (<SDF>) { if (m/^\s*(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(\w+)/) { push @{$mols_mm[$count]}, [$4,$1,$2,$3]; next; } if (m/^\$\$\$\$$/) { unshift @{$mols_mm[$count]}, {}; $count++; } } close SDF; unlink "${name}_confab.sdf"; unlink "$name.mol" if $file_name ne "$name.mol"; #write_molden(@mols_mm, "${name}_confab.xyz"); if (! $no_confab_energy && -x $obminimize_x) { local $method = 'MMFF94'; local $MM_key = '-c 5e-6'; print "\nRun obminimize to get energy of confab's conformers\n"; printf "% 4s%10s\n", 'Number', 'Energy'; my $count = 0; foreach (@mols_mm) { $_ = obminimize2xyz($_); printf "% 04d%10.2f\n", $count++, $_->[0]{Energy}; } print "\n"; } } ########################################################################## # elsif (defined($babel) && $babel ne 'xyz') { # print "Multi molecula file in $babel format\n"; # die "No babel executable $babel_x\n" unless -x $babel_x; # my @args = ($babel_x, "-i$babel", "$name.$ext", "-oxyz", '-m', "babel_$name"); # print "\nRun babel:\n@args\n"; # system(@args) == 0 or do {warn "System\n@args\nfailed: $?\n"; next}; # $ext = 'xyz'; # my @files = glob("babel_$name*"); # @mols_mm = read_molden(@files); # unlink @files; # } # ########################################################################## # elsif ($ext eq 'sdf') { @mols_mm = sdf2mol("$name.$ext"); } elsif ($ext ne 'xyz') { print "$name.$ext: multi molecula file in $ext format\n"; unless (-x $babel_x) { warn "No babel executable $babel_x for conversion $ext to xyz\n"; next; } my @args = ($babel_x, "$name.$ext", "babel_$name.xyz"); print "\nRun babel:\n@args\n"; system(@args) == 0 or do {warn "System\n@args\nfailed: $?\n"; next}; $ext = 'xyz'; next unless -f "babel_$name.xyz"; @mols_mm = read_molden("babel_$name.xyz"); unlink "babel_$name.xyz"; } ########################################################################## else { print "Multi molecula file in xyz format\n"; @mols_mm = read_molden("$name.$ext"); } ########################################################################## ########################################################################## next unless @mols_mm; # Делаем проверку на число атомов. Нумеруем молекулы for (my $i=0; $i<@mols_mm; $i++) { if (@{$mols_mm[$i]} != @{$mols_mm[0]}) { die "Point ".($i+1)." not isomer in $name. Die.\n"; } $mols_mm[$i][0]{'Point'} = $i+1; $mols_mm[$i][0]{'File_name'} = $file_name; } if (@mols_mm > 1) { unless ($no_sort_energy) { print "Sorting..."; @mols_mm = sort {($a->[0]{'Energy'} || 0) <=> ($b->[0]{'Energy'} || 0)} @mols_mm; #@mols_mm = #map {$_->[1]} #sort {($a->[0]) <=> ($b->[0])} #map {[$_->[0]{'Energy'} || 0, $_]} #@mols_mm; } print " done\n"; } if ($only_gen) { print_write_molden(\@mols_mm, $name, undef); exit; } #my $el_f = $ellips ? $ellips : 2; my $el_f = 3; #$el_f = 0 if defined($ellips) && $ellips==0; my $ellips_format = "%.$el_f"."f_%.$el_f"."f_%.$el_f"."f"; my ($small_rms, $all_rms); if (! $no_cf and @mols_mm > 1) { # Вычисляем формат эллипсоида инерции my @eigen_num = tenzor(@mols_mm); for (my $i=0; $i<@mols_mm; $i++) { $mols_mm[$i][0]{'Ellips'} = sprintf $ellips_format, @{$eigen_num[$i]}; } print "\nSorted by energy conformers: ", $#mols_mm+1, "\n"; printf "Number Energy Inertia_ellipsoid\n"; for (my $i=0; $i<@mols_mm; $i++) { printf " %04d %10.4f %.2f_%.2f_%.2f\n", $mols_mm[$i][0]{'Point'}, $mols_mm[$i][0]{'Energy'}||0, @{$eigen_num[$i]}; } if (@mols_mm > 1) { #print "\nУдаление дублей и сортировка:\n"; ($small_rms, $all_rms) = small_RMS(\@mols_mm, $RMS || ''); #print "\nRMS matrix\n"; #foreach (@$all_rms) { # foreach (@$_) { # printf( $_ ? ("%5.2f ",$_) : ("%5s ", '-')); # } # print "\n"; #} #dd $all_rms; print "\nRemoval of doubles and sorting:\n"; foreach (sort {$a<=>$b} keys %$small_rms) { my $I = $1 if $mols_mm[$_][0]{'Point'} =~ /(\d+)/; my $J = $1 if $mols_mm[$small_rms->{$_}[0]][0]{'Point'} =~ /(\d+)/; printf "%04d = %04d, RMS %5.2f\n", $I, $J, $small_rms->{$_}[1]; } print "was ", scalar @mols_mm, ", is "; @mols_mm = delete_same(\@mols_mm,$small_rms); print scalar @mols_mm, " conformers\n"; } } # Проверяем, возможны ли вычисления заданным методом # Если невозможны, печатаем причину и устанавливаем $no_mopac = 1 if (! $no_mopac) { if ($method =~ /RM1|PM6|PM7|QM/) { for (my $i=1; $i<@{$mols_mm[0]}; $i++) { my $element = ucfirst $mols_mm[0][$i][0]; if ($method eq 'RM1' && ! exists $RM1{$element}) { print "No $element in RM1. Force -no_mopac\n"; $no_mopac = 1; last; } elsif ($method eq 'PM6' && ! exists $PM6{$element}) { print "No $element in PM6. Force -no_mopac\n"; $no_mopac = 1; last; } elsif ($method eq 'PM7' && ! exists $PM7{$element}) { print "No $element in PM7. Force -no_mopac\n"; $no_mopac = 1; last; } elsif ($method eq 'QM' && ! exists $QM{$element}) { print "No $element in QM. Force -no_mopac\n"; $no_mopac = 1; last; } } } if ($method =~ /MMFF94s?|UFF|Ghemical|Gaff/ && @mols_mm >1) { unless (obminimize2xyz($mols_mm[0])) { print "No parameters in $method. Force -no_mopac\n"; $no_mopac = 1; last; } } if ($method =~ /tinker/ && @mols_mm >1) { unless (tminimize2xyz($mols_mm[0],"$tinker_params/$tinker_MM",$tinker_rms)) { print "Something wrong for -method=tinker. Force -no_mopac\n"; $no_mopac = 1; last; } } } if (defined $no_mopac) { print_write_molden(\@mols_mm, $name, $all_rms); next; # К следующему файлу } if ($method eq 'QM') { print "\nRun Priroda ($method method):\n$p14_x $name.NUM.in\n"; } elsif ($method =~ /RM1|PM6|PM7|AM1|PM3|MNDO/) { print "\nRun mopac ($method hamiltonian):\n$mopac_x $name.NUM.mop\n"; } elsif ($method =~ /MMFF94s?|UFF|Ghemical|Gaff/) { print "\nRun obminimize ($method force field):\n$obminimize_x -ff $method $MM_key -oxyz $name.NUM.$EXT\n"; } elsif ($method =~ /tinker/) { print "\nRun tinker optimization ($tinker_mm force field):\n"; print "$tminimize_x $name.NUM.txyz $tinker_params/$tinker_MM "; print 'A A ' if $tminimize_x =~ /newton$/; print "$tinker_rms\n"; } printf "%4s %22s%17s%12s\n", 'NUM', 'Point', 'Energy', 'Symmetry'; my @mols_mop; #dd $_->[0] foreach @mols_mm; exit; for (my $count=0; $count<@mols_mm; $count++) { my $mol = $mols_mm[$count]; # Переносим первоначальны?? номер молекулы $mol->[0]{'Point'} = "($mol->[0]{'Point'})" if $mol->[0]{'Point'} =~ /,/; $mols_mop[$count][0]{'Point'} = $mol->[0]{'Point'}; #dd $mols_mop[$count][0]{'Point'}; my $file = $name . '.' . sprintf("%04d", $count+1); if ($calc_charge) { my $sum = 1 - $mult%2; for (my $i=1; $i<@$mol; $i++) { $sum += element($mol->[$i][0]); } $charge = $sum%2; } if ($method eq 'QM') { open IN, '>', "$file.in" or do {warn "Can't write $file.in: $!\n"; next}; print IN "\$system mem=256 disk=-256 \$end\n"; print IN "\$control task=optimize theory=qm_n3 basis=$p14_basis \$end\n"; print IN "\$optimize tol=1e-6 steps=300 cart=1 \$end\n"; print IN "\$molecule charge=$charge mult=$mult cartesian\n"; for (my $i=1; $i<@$mol; $i++) { printf IN " %4d %12.8f %12.8f %12.8f\n", element($mol->[$i][0]), $mol->[$i][1], $mol->[$i][2], $mol->[$i][3]; } print IN "\$end\n"; close IN; if ($keep_out) { open OUT, '>', "$file.out" or die "Can't open : $file.out$!\n"; } my @arg = ($mpiexec_x, '-np', $np, $p14_x, "$file.in"); open L, '-|', @arg or die "Can't run @arg: $!\n"; my $n = @$mol - 1; my $OK; while (<L>) { print OUT if $keep_out; if (m/Atomic Coordinates:/) { for (my $i=1; $i<=$n; $i++) { my $line = <L>; print OUT $line if $keep_out; $mols_mop[$count][$i] = [split ' ', $line]; $mols_mop[$count][$i][0] = element($mols_mop[$count][$i][0]); } } if (m/^\s+Energy:\s+(\S+)/) { $mols_mop[$count][0]{'Energy'} = $1; $mols_mop[$count][0]{'Energy_SP'} ||= $1 if $QM_SP; } if (m/^\s+max. gradient\s+=\s*(\S+)/) { print "\r Gradient $1 Energy $mols_mop[$count][0]{'Energy'}"; } if (m/OPTIMIZATION CONVERGED/) { $OK = 1; } } close L; close L if $keep_out; unlink "$file.in" unless $keep_out; unless ($OK) { printf "\r %04d Failure \n", $count+1; next; } if (-x $symmetry_x) { symmetr($mols_mop[$count]); } else { warn "No $symmetry_x\n"; $mols_mop[$count][0]{'Energy'} = ''; } printf "\r%04d %22s%17.6f%12s\n", $count+1, $mols_mop[$count][0]{'Point'}, $mols_mop[$count][0]{'Energy'} || '', $mols_mop[$count][0]{'Symmetry'} || ''; } elsif ($method =~ /RM1|PM6|PM7|AM1|PM3|MNDO/) { open MOP, '>', "$file.mop" or do {warn "Can't write $file.mop: $!\n"; next}; print MOP " $method CHARGE=$charge $mop_key $mopac_key\n\n\n"; for (my $i=1; $i<@$mol; $i++) { printf MOP " %-2s %12.8f %12.8f %12.8f\n", @{$mol->[$i]}; } close MOP; unlink glob("$file.{arc,out}"); my @args = ($mopac_x, "$file.mop"); #print "\nRun mopac:\n@args\n"; system(@args) == 0 or do {warn "System\n@args\nfailed: $?\n"; next}; if (-e "$file.arc") { open ARC, '<', "$file.arc" or do {warn "Can't open $file.arc: $!\n"; next}; while (<ARC>) { if (m/HEAT OF FORMATION\s+=\s+(\S+)/) { $mols_mop[$count][0]{'Energy'} = $1; } if (m/POINT GROUP:\s+(\S+)/) { $mols_mop[$count][0]{'Symmetry'} = $1; } if (m/FINAL GEOMETRY OBTAINED/) { <ARC>; <ARC>; <ARC>; while (<ARC>) { last if m/^\s*$/; push @{$mols_mop[$count]}, [(split)[0,1,3,5]]; } } } close ARC; # open OUT, '<', "$file.out" or do {warn "Can't open $file.out: $!\n"; next}; # while (<OUT>) { # next if 1../SCF FIELD WAS ACHIEVED/; # if (m/FINAL HEAT OF FORMATION =\s+(\S+)/) { # $mols_mop[$count][0]{'Energy'} = $1; # } # if (m/POINT GROUP:\s+(\S+)/) { # $mols_mop[$count][0]{'Symmetry'} = $1; # } # if (m/CARTESIAN COORDINATES/) { # <OUT>; <OUT>; <OUT>; # while (<OUT>) { # last if m/^\s*$/; # push @{$mols_mop[$count]}, [(split)[1..4]]; # } # } # } # close OUT; rename "$file.arc", "$file.arcSA"; if (! defined $mols_mop[$count][0]{'Energy'}) { printf " %04d Failure\n", $count+1; next; } printf "%04d %22s%17.4f%12s\n", $count+1, $mols_mop[$count][0]{'Point'}, $mols_mop[$count][0]{'Energy'} || '', $mols_mop[$count][0]{'Symmetry'} || ''; } else { warn "No $file.arc. Perhaps, mopac error. Check $file.FAIL.out\n"; #$keep_out = 1; rename("$file.out", "$file.FAIL.out"); rename("$file.mop", "$file.FAIL.mop"); } } elsif ($method =~ /MMFF94s?|UFF|Ghemical|Gaff/) { $mols_mop[$count] = obminimize2xyz($mol); printf "%04d %22s%17.4f%12s\n", $count+1, $mols_mop[$count][0]{'Point'}, $mols_mop[$count][0]{'Energy'} || '', $mols_mop[$count][0]{'Symmetry'} || ''; } elsif ($method =~ /tinker/) { $mols_mop[$count] = tminimize2xyz($mol,"$tinker_params/$tinker_MM",$tinker_rms); printf "%04d %22s%17.4f%12s\n", $count+1, $mols_mop[$count][0]{'Point'}, $mols_mop[$count][0]{'Energy'} || '', $mols_mop[$count][0]{'Symmetry'} || ''; } #write_molden($mols_mop[$count], "$file.xyz"); $mols_mop[$count][0]{'count'} = $count+1; } if ($keep_out) { unlink glob("$name.[0-9][0-9][0-9]*.{den,res}"); } else { unlink grep {! /\.FAIL\./} glob("$name.[0-9][0-9][0-9]*.{mop,out,den,res}"); } @mols_mop = grep {@$_ > 1} @mols_mop; if (@mols_mop > 1) { unless ($no_sort_energy) { @mols_mop = sort {$a->[0]{'Energy'} <=> $b->[0]{'Energy'}} grep {defined $_->[0]{'Energy'}} @mols_mop; } #dd @mols_mop; if (! $no_cf and @mols_mop >1) { my @eigen_num = tenzor(@mols_mop); for (my $i=0; $i<@mols_mop; $i++) { $mols_mop[$i][0]{'Ellips'} = sprintf $ellips_format, @{$eigen_num[$i]}; } (my $small_rms, $all_rms) = small_RMS(\@mols_mop, $RMS || ''); print "\nRemoval of doubles and sorting:\n"; #dd $_->[0] foreach @mols_mop; exit; foreach (sort {$a<=>$b} keys %$small_rms) { my $I = $1 if $mols_mop[$_][0]{'Point'} =~ /(\d+)/; my $J = $1 if $mols_mop[$small_rms->{$_}[0]][0]{'Point'} =~ /(\d+)/; printf "%04d = %04d, RMS %5.2f\n", $I, $J, $small_rms->{$_}[1]; } # foreach (sort {$a<=>$b} keys %$small_rms) { # printf "%04d = %04d, RMS %5.2f\n", # $mols_mop[$_][0]{'Point'}, # $mols_mop[$small_rms->{$_}[0]][0]{'Point'}, ## $_+1, ## $small_rms->{$_}[0]+1, # $small_rms->{$_}[1]; # } print "was ", scalar @mols_mop, ", is "; @mols_mop = delete_same(\@mols_mop,$small_rms); print scalar @mols_mop, " conformers\n"; } } print_write_molden(\@mols_mop, $name, $all_rms) if @mols_mop > 0; for (my $i=0; $i<@mols_mop; $i++) { rename sprintf("$name.%04d.arcSA",$mols_mop[$i][0]{'count'}), sprintf("$name.%04d.arc",$i+1); } unlink glob("$name.*.arcSA"); close LOG if $logfile; } sub print_write_molden { my @mols = @{$_[0]}; my $name = $_[1]; my $all_rms = $_[2]; if (@mols==1) { print "\nWrite $name.0001.xyz\n"; write_molden($mols[0], "$name.0001.xyz"); return; } print "\nWrite $name.NNNN.xyz\n"; my $maxL = 0; foreach (@mols) { my $l = length $_->[0]{'Point'}; $maxL = $l if $maxL < $l; } my $k = 1; printf "%${maxL}s %3s%12s%10s%20s\n", '', 'NNNN', 'Energy', 'Symmetry', 'Inertia_ellipsoid'; foreach my $mol (@mols) { my $NNNN = sprintf("%04d",$k++); printf $mol->[0]{'Point'} ? ("%${maxL}s => ", $mol->[0]{'Point'}) : ("%${maxL}s", ' '); printf "%04d", $NNNN; printf "%12.4f", $mol->[0]{'Energy'}||0; printf "%10s", $mol->[0]{'Symmetry'}||'-'; printf "%20s", $mol->[0]{'Ellips'}||'-'; print "\n"; write_molden($mol, "$name.$NNNN.xyz") } my @similar; if ($all_rms) { print "\nRMS matrix in (RMS = $RMS) units\n"; print 'NNNN|', (map {$_%10} 1..$#mols+1), "\n"; #print "---|", '-'x@mols, "\n"; for (my $i=0; $i<@mols; $i++) { printf "%04d|", $i+1; for (my $j=0; $j<@mols; $j++) { print(' '),next if $i==$j; #print(0),next if $i==$j; my $I = $1 if $mols[$i][0]{'Point'} =~ /(\d+)/; my $J = $1 if $mols[$j][0]{'Point'} =~ /(\d+)/; my $rms = $all_rms->[$I][$J] || $all_rms->[$J][$I]; if ($rms =~ /\d+(\.\d+)?/) { $rms = int($rms/$RMS); if ($rms == 1 && $i < $j) { push @similar, [sprintf("%04d", $i+1), sprintf("%04d", $j+1)]; } $rms = 9 if $rms > 9; } print $rms; } print "\n"; } } if (@similar) { print "\nSimilar structures:\n"; foreach (@similar) { my @a = map {"$name.$_.xyz"} @$_; print "@a\n"; } } } sub small_RMS { ####################################################################### # Первый аргумент -- ссылка на массив молекул # следом м.б. минимальный RMS # Возвращается хеш ??ндексов молекул, которые следует удалить из массива ####################################################################### my @mols = @{$_[0]}; my $rms_limit = $_[1] || 0.1; my (%small_rms, @all_rms); # return my $diagonal = $#mols; $diagonal = $#mols if $diagonal > $#mols; print "\nRMSD matrix by diagonals in (RMS = $rms_limit) units\n"; for (my $d=0; $d<$diagonal; $d++) { my $count=1; for (my $i=0; $i<$#mols-$d; $i++) { my $ai = $1 if $mols[$i][0]{'Point'} =~ /(\d+)/; my $aj = $1 if $mols[$i+$d+1][0]{'Point'} =~ /(\d+)/; if (exists($small_rms{$i+$d+1}) || exists($small_rms{$i})) { #$all_rms[$i][$i+$d] = '-'; $all_rms[$ai][$aj] = '-'; print '-'; next; } $count++; if ($energy && defined($mols[$i][0]{'Energy'}) && defined($mols[$i+$d+1][0]{'Energy'}) && abs($mols[$i][0]{'Energy'}-$mols[$i+$d+1][0]{'Energy'}) > $energy) { #$all_rms[$i][$i+$d] = 'E'; $all_rms[$ai][$aj] = 'E'; print 'E'; next; } if (defined($ellips) && defined($mols[$i][0]{'Ellips'}) && defined($mols[$i+$d+1][0]{'Ellips'}) && rms_ellips($mols[$i][0]{'Ellips'}, $mols[$i+$d+1][0]{'Ellips'})>$ellips) { $all_rms[$ai][$aj] = 'L'; print 'L'; next; } my $rms = tenzor_sort(@mols[$i,$i+$d+1], $sort_RMS); #printf STDERR "%3d: %3d vs %-3d", $d+1, $i+1,$i+$d+2; my $print_rms = int($rms/$rms_limit); $print_rms = 9 if $print_rms > 9; print $print_rms; #$all_rms[$i][$i+$d] = $rms; $all_rms[$ai][$aj] = $rms; if ($rms < $rms_limit) { $small_rms{$i+$d+1} = [$i,$rms]; #$mols[$i+$d+1][0]{'Point'} .= ','.$mols[$i][0]{'Point'}; $mols[$i][0]{'Point'} .= ','.$mols[$i+$d+1][0]{'Point'}; #print STDERR ": equal" } #print STDERR "\r", ' 'x25, "\r"; #print STDERR "\n"; } #printf STDERR "\r% 7d\r", $diagonal-$d+1; print "\n"; } return \%small_rms, \@all_rms; } sub delete_same { ####################################################################### # Первый аргумент -- ссылка на массив молекул # Второй -- ссылка на двумерный массив RMS (следом м.б. минимальный RMS) # либо ссылка на хеш индексов молекул, которые следует удалить из массива. # Возвращается отсортированный по энергии массив молекул без дублей ####################################################################### my @mols = @{$_[0]}; #print 'Was ', $#mols+1; if (ref($_[1]) eq 'ARRAY') { my @rms = @{$_[1]}; my $rms_limit = $_[2] || 0.1; # for (my $i=0; $i<@rms; $i++) { # for (my $j=0; $j<@{$rms[$i]}; $j++) { # next unless defined $rms[$i][$j]; # if ($rms[$i][$j] < $rms_limit) { # $mols[$i][0]{'Point'} .= ','.$mols[$j][0]{'Point'}; # $mols[$j][0]{'Point'} .= ','.$mols[$i][0]{'Point'}; # } # } # } my $is_null_rms = sub { my $i = shift; my @rms = @_; for (my $j=0; $j<@rms; $j++) { next unless defined $rms[$i][$j]; return 0 if $rms[$i][$j] < $rms_limit; } return 1 }; my $k = 0; @mols = sort {($a->[0]{'Energy'} || 0) <=> ($b->[0]{'Energy'} || 0)} @mols unless $no_sort_energy; @mols = grep {$is_null_rms->($k++,@rms)} @mols; } elsif (ref($_[1]) eq 'HASH') { my %small_rms = %{$_[1]}; my $k = 0; @mols = sort {($a->[0]{'Energy'} || 0) <=> ($b->[0]{'Energy'} || 0)} @mols unless $no_sort_energy; @mols = grep {! exists $small_rms{$k++}} @mols; } #print '; Is ', $#mols+1, " conformers\n"; return @mols; } sub tenzor { ########################################################################### ## Принимает массив ссылок (молекул) и модифицирует их по месту, ## перемещая в це??тр масс и ориентируя по осям тензора инерции. ## Использует внешний хэш %massa и функцию jacobi(). ########################################################################### my @eigen_num; print "\nOrientation of molecula by inertia tenzor\n"; foreach my $mol (@_) { my (@m,@atom,@x,@y,@z); my $N = $#{$mol}; # Создаем массив атомных масс for (my $i=1; $i<=$N; $i++) { $atom[$i] = $mol->[$i][0]; $x[$i] = $mol->[$i][1]; $y[$i] = $mol->[$i][2]; $z[$i] = $mol->[$i][3]; $m[$i] = $massa{$atom[$i]}; } # Вычисляем молекулярную массу my $M = sum(@m[1..$N]); #print "$M\n"; # Перемещаем геометрию в центр масс my $sum; #&write_molden($mol, 'ttt'); $sum = sum(map { $x[$_]*$m[$_] } 1..$N)/$M; $_ -= $sum for @x[1..$N]; $sum = sum(map { $y[$_]*$m[$_] } 1..$N)/$M; $_ -= $sum for @y[1..$N]; $sum = sum(map { $z[$_]*$m[$_] } 1..$N)/$M; $_ -= $sum for @z[1..$N]; # Суть: молекулу представляем в виде эллипсоида инерции. # Длина главных осей этого эллипсоида определяются собственными значениями # тензора инерции, а их направление - соответствующими собственными векторами. # Г. Голдстейн. Классическая механика. Глава 5. # Вычисляем тензор инерции # | Ixx Ixy Ixz | # | Ixy Iyy Iyz | # | Ixz Iyz Izz | my (@r2,$Ixx,$Ixy,$Ixz,$Iyy,$Iyz,$Izz); $r2[$_] = $x[$_]**2+$y[$_]**2+$z[$_]**2 for 1..$N; $Ixx = sum(map {$m[$_]*($r2[$_]-$x[$_]**2)} 1..$N); $Iyy = sum(map {$m[$_]*($r2[$_]-$y[$_]**2)} 1..$N); $Izz = sum(map {$m[$_]*($r2[$_]-$z[$_]**2)} 1..$N); $Ixy = -sum(map {$m[$_]*$x[$_]*$y[$_]} 1..$N); $Ixz = -sum(map {$m[$_]*$x[$_]*$z[$_]} 1..$N); $Iyz = -sum(map {$m[$_]*$y[$_]*$z[$_]} 1..$N); #($Ixx,$Ixy,$Ixz,$Iyy,$Iyz,$Izz) = (6,-2,2,5,0,7); #print "Тензор инерции:\n"; #printf "%12.8f %12.8f %12.8f\n%12.8f %12.8f %12.8f\n%12.8f %12.8f %12.8f\n", # $Ixx, $Ixy, $Ixz, $Ixy, $Iyy, $Iyz, $Ixz, $Iyz, $Izz; my @eigen = jacobi($Ixx,$Ixy,$Ixz,$Iyy,$Iyz,$Izz) or die "Eigen error\n"; push @eigen_num, [sqrt($eigen[0][0]/$M), sqrt($eigen[1][0]/$M), sqrt($eigen[2][0]/$M)]; # for debug # for (my $i=0; $i<=2; $i++) { # printf STDERR "%10.6f: %10.4f %10.4f %10.4f\n", $eigen_num[-1][$i], @{$eigen[$i][1]}; # } # print STDERR "\n"; # Преобразуем геометрию к новой системе координат (собственные вектора) my @new_geom; for (my $i=0; $i<$N; $i++) { for (0..2) { $new_geom[$_][$i] = $x[$i+1]*$eigen[$_][1][0] + $y[$i+1]*$eigen[$_][1][1] + $z[$i+1]*$eigen[$_][1][2]; } } #write_molden(); # Переставляем оси эллипсоида, чтобы хорошо смотрелось в молдене @x[1..$N] = @{$new_geom[1]}; @y[1..$N] = @{$new_geom[2]}; @z[1..$N] = @{$new_geom[0]}; # Из новой геометрии делаем две молекулы (одну с отражением по z) my ($mol1,$mol2); for (my $i=1; $i<=$N; $i++) { $mol2->[$i][0] = $mol1->[$i][0] = $mol->[$i][0]; $mol2->[$i][1] = $mol1->[$i][1] = $x[$i]; $mol2->[$i][2] = $mol1->[$i][2] = $y[$i]; $mol1->[$i][3] = $z[$i]; $mol2->[$i][3] = -$z[$i]; } my $znakz = 1; #my $rms1 = get_rms($mol,$mol1); my $rms2 = get_rms($mol,$mol2); print '='; #printf "%20.10f %20.10f\n", $rms1, $rms2; $znakz = -1 if $rms2 < 1e-8; # Перезаписываем молекулу for (my $i=1; $i<=$N; $i++) { $mol->[$i][1] = $x[$i]; $mol->[$i][2] = $y[$i]; $mol->[$i][3] = $znakz*$z[$i]; } #dd $mol; # Проверка на вырожденность собственных чисел тензора инерции # ... по плоскости xy #dd @eigen_num; if ($eigen_num[-1][1]-$eigen_num[-1][2] < 0.001) { #dd $mol; my @d = grep {$_->[0] > 0.1} proections($mol,'xy'); #dd @d; # Если есть хотя бы 3 равноудаленных атома... if (@d>2 && $d[1][0]-$d[0][0]<0.001 && $d[2][0]-$d[0][0]<0.001) { # будем переориентировать, направив ось x к первому из них my $n = $d[0][1]; my $cos = $mol->[$n][1]/sqrt($mol->[$n][1]**2+$mol->[$n][2]**2); my $sin = sqrt(1-$cos**2); $sin = -$sin if $mol->[$n][2] < 0; for (my $i=1; $i<@$mol; $i++) { my $x = $cos*$mol->[$i][1] + $sin*$mol->[$i][2]; my $y = -$sin*$mol->[$i][1] + $cos*$mol->[$i][2]; $mol->[$i][1] = $x; $mol->[$i][2] = $y; } } } # ... по плоскости xz if ($eigen_num[-1][0]-$eigen_num[-1][1] < 0.001) { #dd $mol; my @d = grep {$_->[0] > 0.1} proections($mol,'xz'); #dd @d; # Если есть хотя бы 3 равноудаленных атома... if (@d>2 && $d[1][0]-$d[0][0]<0.001 && $d[2][0]-$d[0][0]<0.001) { # будем переориентировать, направив ось x к первому из них my $n = $d[0][1]; my $cos = $mol->[$n][1]/sqrt($mol->[$n][1]**2+$mol->[$n][3]**2); my $sin = sqrt(1-$cos**2); $sin = -$sin if $mol->[$n][2] < 0; for (my $i=1; $i<@$mol; $i++) { my $x = $cos*$mol->[$i][1] + $sin*$mol->[$i][3]; my $z = -$sin*$mol->[$i][1] + $cos*$mol->[$i][3]; $mol->[$i][1] = $x; $mol->[$i][3] = $z; } } #dd $mol; } } print "\n"; return @eigen_num; } sub tenzor_sort { ########################################################################### ## Принимает две ссылки (молекулы) и sort_RMS третьим параметром. ## Молекулы должны быть предварительно пропущены через tenzor. ## Возвращает rms между этими молекулами. ## Если rms < sort_RMS, то порядок атомов второй молекулы сортируется по первой. ## rms = sqrt(sum(m*r**2)/sum(m)) ## Использует функции znak() и kinen() ########################################################################### my ($mol0, $mol, $sort_RMS) = @_; $sort_RMS ||= 0; my $N = $#{$mol0}; # Создаем структуру данных "молекула", # объединяющую символы атомов, их порядок в геометрии и координаты. # Это хэш (атомных символов) хэшей (их номера от 1 до N) массивов (координат) my (%mol0, %mol); foreach (1..$N) { my @a = @{$mol0->[$_]}; $mol0{$a[0]}->{$_} = [@a[1..$#a]]; @a = @{$mol->[$_]}; $mol{$a[0]}->{$_} = [@a[1..$#a]]; } # Нужно согласовать направления осей текущей # геометрии с направлениями предыдущей. Эти направления после уже # проделанных преобразования могут отличаться только знаками. Сложность # же заключается в том, что порядок нумерации атомов в двух геометриях # может быть различен. Но предполагаем, что две соседние геометрии # не сильно отличаются друг от друга # Вычисляем "кинетич. энергию" между текущей и предыдущей молекулами # В возвращаемом массиве 0-й элемент - это "кинетич. энергия", остальные - # пермутационный вектор, оп??еделяющий порядок нумерации текущей молекулы, # при котором она лучше всего накладывается на предыдущую. my @jmin = kinen(\%mol,\%mol0); # Далее будем последовательно менять знаки координат текущей молекулы # (отражение) и смотреть, при каком отражении будет мин. кин. энергия # Последовательность отражений: # 0 +x +y +z # 1 +x -y +z # 2 -x -y +z # 3 -x +y +z # 4 -x +y -z # 5 +x +y -z # 6 +x -y -z # 7 -x -y -z my $reflect = 0; # Счетчик отражений my $reflectmin = 0; # Номер отражения с минимальн. кинетич. энергией # Для последовательности отражений, подобранной так, чтобы меняя каждый раз # знак только одной координаты (0 - x, 1 - y, 2 - z), охватить все варианты my @j; foreach (1,0,1,2,0,1,0) { # Меняем знак соответствующей координаты znak(\%mol, $_); $reflect++; # Счетчик отражений # Если нужно сохранить конфигурацию, пропускаем вычисления на нечетных шагах # (останутся только повороты вокруг осей z,x,y) # 0 +x +y +z # 2 -x -y +z # 4 -x +y -z # 6 +x -y -z unless ($no_ret) { next if $reflect%2; } # Вычисляем кин. энергию @j = kinen(\%mol,\%mol0); #print "@j\n"; # Смотрим, не лучше ли очередная кин. энергия if ($j[0] < $jmin[0]) { @jmin = @j; $reflectmin = $reflect; } } #warn "$reflectmin @jmin\n"; # Вспомогательный ма??сив, определяющий, как нужно поменять знаки координат, # чтобы из молекулы, полученной после последнего отражения (-x -y -z), # восстановить молекулу после i-го отражения my @ar = ([0,1,2],[0,2],[2],[1,2],[1],[0,1],[0],[]); # Меняем знаки, чтобы получить наименьшую кин. энергию znak(\%mol, @{$ar[$reflectmin]}); # Переносим полученные координаты из "молекулы" %mol в массивы @x,@y,@z # в зависимости от значения опции -sort if ($sort_RMS > $jmin[0]) { # Сортировка по предыдущей геометрии my $i = 1; #print "@jmin\n"; foreach my $j (@jmin[1..$#jmin]) { foreach my $atom (keys %mol) { if (exists $mol{$atom}{$j}) { @{$mol->[$i]} = ($atom, @{$mol{$atom}{$j}}); $i++; } } } } return $jmin[0]; } sub znak { ########################################################################### ## Менять знаки координат в структуре $h{$atom}{$index}[$xyz] ## Usage: znak(\%h, 0|1|2); x - 0, y - 1, z - 2 ########################################################################### my %h = %{$_[0]} and shift; foreach my $xyz (@_) { foreach my $atom (keys %h) { for (keys %{$h{$atom}}) { $h{$atom}{$_}[$xyz] *= -1; } } } } sub kinen { ########################################################################### ## Вычисляет "кинетич. энергию" между сравниваемой и п??едыдущей молекулами ## Принимает две ссылки на сравниваемую и предыдущую молекулы. Молекула - ## это хэш (атомных символов) хэшей (их номера от 1 до N) массивов (координат) ## В возвращаемом массиве 0-й элемент - это "кинетич. энергия", остальные - ## пермутационный вектор, определяющий порядок нумерации текущей молекулы, ## при котором она лучше всего накладывается на предыдущую. ## "кинетич. энергия" rms = sqrt(sum(m*r**2)/sum(m)) ## Использует внешний хэш %massa ########################################################################### # Делаем независимую копию 1-й (сравниваемой) молекулы my %mol; foreach my $atom (keys %{$_[0]}) { while (my ($ind,$aref) = each %{$_[0]{$atom}}) { $mol{$atom}{$ind} = [@{$aref}]; } } # 2-я (предыдущая) молекула my %mol0 = %{$_[1]}; # Вычисляем молекулярную массу my $M = sum(map {$massa{$_}*(keys %{$mol{$_}})} keys %mol); my $E = 0; my @j; # Для углеродов, водородов и т.д. начиная с более тяжелых foreach my $atom (sort {$massa{$b}<=>$massa{$a}} keys %mol0) { # Если в молекулах разное кол-во, например, углеродов, то выход if (scalar(keys %{$mol{$atom}}) != scalar(keys %{$mol0{$atom}})) { warn "Not isomers ($atom)\n"; return; } # Для каждого из атомов данного сорта 2-й молекулы # отсортированных по увеличению расстояния до центра масс foreach my $i (map {$_->[0]} sort {$a->[1]<=>$b->[1]} map {[$_,$mol0{$atom}{$_}[0]**2+$mol0{$atom}{$_}[1]**2+$mol0{$atom}{$_}[2]**2]} keys %{$mol0{$atom}}) { my $rmin = 1e10; my $jmin; # перебираем атомы того же сорта 1-й молекулы foreach my $j (keys %{$mol{$atom}}) { # вычисляем расстояние между атомами двух молекул #my $r = sum(map {($mol{$atom}{$j}[$_] - $mol0{$atom}{$i}[$_])**2} 0..2); my $r = ($mol{$atom}{$j}[0] - $mol0{$atom}{$i}[0])**2 + ($mol{$atom}{$j}[1] - $mol0{$atom}{$i}[1])**2 + ($mol{$atom}{$j}[2] - $mol0{$atom}{$i}[2])**2; # запоминая ближайший атом 1-й молекулы if ($r<$rmin) { $rmin = $r; $jmin = $j; } } # Сумимируем кинетическую энергию $E += $rmin*$massa{$atom}; # i-му атому 2-й молекулы ставим в соответствие ближайший атом из 1-й $j[$i] = $jmin; # и удаляем его из 1-й молекулы delete $mol{$atom}{$jmin}; } # Записываем накопленную кин. энергию в начало массива $j[0] = $E; } $j[0] = sqrt($j[0]/$M); # Возвращаем массив return @j; } sub jacobi { ########################################################################### ## Собственные числа и векторы симметричной матрицы методом Якоби. ## Код из Жориной программы dte (C). ## Принимает массив ($Ixx,$Ixy,$Ixz,$Iyy,$Iyz,$Izz). ## Возвращает отсортированный (по убыванию собств. чисел) массив двухэлементных ## массивов. 1-й элемент - собств. число, 2-й - соответствующий ему вектор. ########################################################################### my @a = ([ $_[0],$_[1],$_[2] ], [ $_[1],$_[3],$_[4] ], [ $_[2],$_[4],$_[5] ]); my ($j, $iq, $ip, $i, $n); my ($tresh, $theta, $tau, $t, $sm, $s, $h, $g, $c); my $maxit = 50; $n = 3; my @v = ([1,0,0], [0,1,0], [0,0,1]); my @b = my @d = ($a[0][0],$a[1][1],$a[2][2]); my @z = (0,0,0); my $nrot = 0; for ($i=0; $i<$maxit; $i++) { $sm = abs($a[0][1]) + abs($a[0][2]) + abs($a[1][2]); if ($sm == 0) { return sort {$b->[0] <=> $a->[0]} ([$d[0], [$v[0][0],$v[1][0],$v[2][0]]], [$d[1], [$v[0][1],$v[1][1],$v[2][1]]], [$d[2], [$v[0][2],$v[1][2],$v[2][2]]]); } # The normal return, which relies on # quadratic convergence to machine underflow. $tresh = ($i < 4) ? 0.2*$sm/($n*$n) : 0; # ...on the first three sweeps. # ...thereafter. for ($ip=0; $ip<$n-1; $ip++) { for ($iq=$ip+1; $iq<$n; $iq++) { $g = 100*abs($a[$ip][$iq]); # After four sweeps, skip the rotation if the off-diagonal element is small. if ($i>4 && abs($d[$ip])+$g == abs($d[$ip]) && abs($d[$iq])+$g == abs($d[$iq])) { $a[$ip][$iq] = 0; } elsif (abs ($a[$ip][$iq]) > $tresh) { $h = $d[$iq]-$d[$ip]; if (abs($h)+$g == abs($h)) { $t = $a[$ip][$iq]/$h; } # $t = 1/(2theta) else { $theta = 0.5*$h/$a[$ip][$iq]; $t = 1/(abs($theta)+sqrt(1+$theta*$theta)); if ($theta < 0) {$t = -$t;} } $c = 1/sqrt(1+$t*$t); $s = $t*$c; $tau = $s/(1+$c); $h = $t*$a[$ip][$iq]; $z[$ip] -= $h; $z[$iq] += $h; $d[$ip] -= $h; $d[$iq] += $h; $a[$ip][$iq] = 0; for ($j=0; $j<$ip; $j++) { # Case of rotations 0 <= $j < p. $g = $a[$j][$ip]; $h = $a[$j][$iq]; $a[$j][$ip] = $g-$s*($h+$g*$tau); $a[$j][$iq] = $h+$s*($g-$h*$tau); } for ($j=$ip+1; $j<$iq; $j++) { # Case of rotations p < $j < q. $g = $a[$ip][$j]; $h = $a[$j][$iq]; $a[$ip][$j] = $g-$s*($h+$g*$tau); $a[$j][$iq] = $h+$s*($g-$h*$tau); } for ($j=$iq+1; $j<$n; $j++) { # Case of rotations q < $j < $n. $g = $a[$ip][$j]; $h = $a[$iq][$j]; $a[$ip][$j] = $g-$s*($h+$g*$tau); $a[$iq][$j] = $h+$s*($g-$h*$tau); } for ($j=0; $j<$n; $j++) { $g = $v[$j][$ip]; $h = $v[$j][$iq]; $v[$j][$ip] = $g-$s*($h+$g*$tau); $v[$j][$iq] = $h+$s*($g-$h*$tau); } ++$nrot; } } } for ($ip=0; $ip<$n; $ip++) { $b[$ip] += $z[$ip]; $d[$ip] = $b[$ip]; # Update $d with the sum of tapq, $z[$ip] = 0; # and reinitialize $z. } } return undef; } sub read_molden { ############################################################################ ## Молекула -- ссылка на массив, 0-й элемент -- свойства, следующие - атомы. ## Свойства -- ссылка на хэш с ключами Energy, Symmetry ## Атом -- ссылка на массив [atom, x, y, z, ppm] (ppm может не быть). ## ## Читает xyz. Параметры - имена xyz-файлов. Если параметров нет, то <>. ## Возвращает массив найденных молекул. ############################################################################ local @ARGV = @_ ? @_ : @ARGV; my $num = qr/-?\d+(?:\.\d+)?(?:[eE][+-]?\d+)?/; my @mols; my $line; LOOP: while ($line || defined($line = <>)) { #print $line; if ($line =~/^\s*(\d+)\s*$/) { my @mol; my $N = $1; last LOOP if eof(); next LOOP if eof(ARGV); $line = <>; ($mol[0]{'Energy'}) = $line =~ /(?:^|)\s($num)(?:\s|$)/o; ($mol[0]{'Symmetry'}) = $line =~ /symm\S*\s+(\S+)/i; for (my $i=1; $i<=$N; $i++) { last LOOP if eof(); next LOOP if eof(ARGV); $line = <>; #print $line; if ($line =~ /^\s*([A-Z]{1,2})\s+($num)\s+($num)\s+($num)\s*(.*)/io) { $mol[$i] = [$1,$2,$3,$4,$5]; } else { next LOOP; } } push @mols, \@mol; last LOOP if eof(); } else { undef $line; } } return @mols; } sub write_molden { my $file = pop @_; open F, '>', $file or do {warn "Can't write to $file: $!\n"; return}; foreach my $mol (@_) {; my $N = $#{$mol}; print F " $N\n"; print F " Energy $mol->[0]{'Energy'}" if $mol->[0]{'Energy'}; print F " Energy_SP $mol->[0]{'Energy_SP'}" if $mol->[0]{'Energy_SP'}; print F " Charge $charge" if $charge && $charge != 0; print F " Mult $mult" if $mult && $mult != 1; print F " Point $mol->[0]{'Point'}" if $mol->[0]{'Point'}; print F " Symmetry $mol->[0]{'Symmetry'}" if $mol->[0]{'Symmetry'}; print F " Ellips $mol->[0]{'Ellips'}" if $mol->[0]{'Ellips'}; print F "\n"; for (my $i=1; $i<=$N; $i++) { printf F " %-2s %12.8f %12.8f %12.8f\n", @{$mol->[$i]}; } } close F; } sub element { ############################################################################ ## Возвращает атомный номер по символу либо символ по номеру ## или undef при отсутствии элемента в таблице Менделеева ############################################################################ # Таблица Менделеева 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 $atom = shift; if ($atom=~/^\d+$/) { return undef if $atom > $#ATOM; return $ATOM[$atom]; } $atom = ucfirst $atom; my %ATOM; @ATOM{@ATOM} = 0..$#ATOM; return $ATOM{$atom} if exists $ATOM{$atom}; return undef; } sub sum { my $sum = 0; foreach (@_) { $sum += $_; } return $sum; } sub get_rms { ############################################################################ ## Переделанная на перл coord.c Г.Сальникова ## Принимает две ссылки на геометрии (см. read_molden() из conformers) ## Возвращает RMSD между геометриями ## Во вторую ссылку записывает новую, сглаженную геометрию ## Использует глобальный хэш %massa с массами элементов (см. smooth) ############################################################################ my ($mol0, $mol1) = @_; my ($i, $n, $rot); my (@m, @x0, @y0, @z0, @x1, @y1, @z1, $xc, $yc, $zc, $e, $e0, $mtot, $tg1, $tg2, $phi, $phix, $phiy, $phiz); my $M_PI = atan2(0,-1); my $FLT_EPSILON = 1e-8; my $DEBUG = 0; #/*** Input data preparation ***/ $n = $#{$mol0}; for ($i=0; $i<$n; $i++) { $m[$i] = $massa{$mol0->[$i+1][0]}; die "Inconsistent data\n" if $massa{$mol1->[$i+1][0]} != $m[$i]; $mtot += $m[$i]; $x0[$i] = $mol0->[$i+1][1]; $y0[$i] = $mol0->[$i+1][2]; $z0[$i] = $mol0->[$i+1][3]; $x1[$i] = $mol1->[$i+1][1]; $y1[$i] = $mol1->[$i+1][2]; $z1[$i] = $mol1->[$i+1][3]; } #/*** 1-st translation to center of mass ***/ $xc = $yc = $zc = 0; for ($i=0; $i<$n; $i++) { $xc += $m[$i]*$x0[$i]; $yc += $m[$i]*$y0[$i]; $zc += $m[$i]*$z0[$i]; } $xc /= $mtot; $yc /= $mtot; $zc /= $mtot; for ($i=0; $i<$n; $i++) { $x0[$i] -= $xc; $y0[$i] -= $yc; $z0[$i] -= $zc; } if ($DEBUG) { printf ("1-st molecule in center of mass\n"); for ($i=0; $i<$n; $i++) { printf "%2.0f %10f %10f %10f\n", $m[$i], $x0[$i], $y0[$i], $z0[$i]; } printf "1-st center of mass translation on: (%10f, %10f, %10f)\n", -$xc, -$yc, -$zc; } #/*** 2-nd translation to center of mass ***/ $xc = $yc = $zc = 0; for ($i=0; $i<$n; $i++) { $xc += $m[$i]*$x1[$i]; $yc += $m[$i]*$y1[$i]; $zc += $m[$i]*$z1[$i]; } $xc /= $mtot; $yc /= $mtot; $zc /= $mtot; for ($i=0; $i<$n; $i++) { $x1[$i] -= $xc; $y1[$i] -= $yc; $z1[$i] -= $zc; } if ($DEBUG) { printf ("2-nd molecule in center of mass\n"); for ($i=0; $i<$n; $i++) { printf "%2.0f %10f %10f %10f\n", $m[$i], $x1[$i], $y1[$i], $z1[$i]; } printf "2-nd center of mass translation on: (%10f, %10f, %10f)\n", -$xc, -$yc, -$zc; } $e = 0; for ($i=0; $i<$n; $i++) { $e += $m[$i]*(($x1[$i]-$x0[$i])*($x1[$i]-$x0[$i])+ ($y1[$i]-$y0[$i])*($y1[$i]-$y0[$i])+ ($z1[$i]-$z0[$i])*($z1[$i]-$z0[$i])); } printf ("2*energy: %g\n", $e) if $DEBUG; for ($rot=1; ; $rot++) { $e0 = $e; #/*** Rotation around X ***/ $tg1 = $tg2 = 0; for ($i=0; $i<$n; $i++) { $tg1 += $m[$i]*($y0[$i]*$z1[$i]-$z0[$i]*$y1[$i]); $tg2 += $m[$i]*($y0[$i]*$y1[$i]+$z0[$i]*$z1[$i]); } $phi = atan2 ($tg1, $tg2); for ($i=0; $i<$n; $i++) { $yc = $y1[$i]*cos($phi)+$z1[$i]*sin($phi); $zc = $z1[$i]*cos($phi)-$y1[$i]*sin($phi); $y1[$i] = $yc; $z1[$i] = $zc; } $e = 0; for ($i=0; $i<$n; $i++) { $e += $m[$i]*(($x1[$i]-$x0[$i])*($x1[$i]-$x0[$i])+ ($y1[$i]-$y0[$i])*($y1[$i]-$y0[$i])+ ($z1[$i]-$z0[$i])*($z1[$i]-$z0[$i])); } if ($DEBUG) { for ($i=0; $i<$n; $i++) { printf "%2.0f %10f %10f %10f\n", $m[$i], $x1[$i], $y1[$i], $z1[$i]; } printf "after %d-th rotation around X on: %g (%g deg)\n2*energy: %g\n", $rot, $phi, $phi*180/$M_PI, $e; } $phix = $phi; #/*** Rotation around Y ***/ $tg1 = $tg2 = 0; for ($i=0; $i<$n; $i++) { $tg1 += $m[$i]*($z0[$i]*$x1[$i]-$x0[$i]*$z1[$i]); $tg2 += $m[$i]*($z0[$i]*$z1[$i]+$x0[$i]*$x1[$i]); } $phi = atan2 ($tg1, $tg2); for ($i=0; $i<$n; $i++) { $xc = $x1[$i]*cos($phi)-$z1[$i]*sin($phi); $zc = $z1[$i]*cos($phi)+$x1[$i]*sin($phi); $x1[$i] = $xc; $z1[$i] = $zc; } $e = 0; for ($i=0; $i<$n; $i++) { $e += $m[$i]*(($x1[$i]-$x0[$i])*($x1[$i]-$x0[$i])+ ($y1[$i]-$y0[$i])*($y1[$i]-$y0[$i])+ ($z1[$i]-$z0[$i])*($z1[$i]-$z0[$i])); } if ($DEBUG) { for ($i=0; $i<$n; $i++) { printf "%2.0f %10f %10f %10f\n", $m[$i], $x1[$i], $y1[$i], $z1[$i]; } printf "after %d-th rotation around Y on: %g (%g deg)\n2*energy: %g\n", $rot, $phi, $phi*180/$M_PI, $e; } $phiy = $phi; #/*** Rotation around Z ***/ $tg1 = $tg2 = 0; for ($i=0; $i<$n; $i++) { $tg1 += $m[$i]*($x0[$i]*$y1[$i]-$y0[$i]*$x1[$i]); $tg2 += $m[$i]*($x0[$i]*$x1[$i]+$y0[$i]*$y1[$i]); } $phi = atan2 ($tg1, $tg2); for ($i=0; $i<$n; $i++) { $xc = $x1[$i]*cos($phi)+$y1[$i]*sin($phi); $yc = $y1[$i]*cos($phi)-$x1[$i]*sin($phi); $x1[$i] = $xc; $y1[$i] = $yc; } $e = 0; for ($i=0; $i<$n; $i++) { $e += $m[$i]*(($x1[$i]-$x0[$i])*($x1[$i]-$x0[$i])+ ($y1[$i]-$y0[$i])*($y1[$i]-$y0[$i])+ ($z1[$i]-$z0[$i])*($z1[$i]-$z0[$i])); } if ($DEBUG) { for ($i=0; $i<$n; $i++) { printf "%2.0f %10f %10f %10f\n", $m[$i], $x1[$i], $y1[$i], $z1[$i]; } printf "after %d-th rotation around Z on: %g (%g deg)\n2*energy: %g\n", $rot, $phi, $phi*180/$M_PI, $e; } $phiz = $phi; last if $rot > 100; last if (abs($phix) < $FLT_EPSILON && abs($phiy) < $FLT_EPSILON && abs($phiz) < $FLT_EPSILON && ($e == $e0 || abs($e+$e0) < $FLT_EPSILON || (abs(($e-$e0)/($e+$e0)) < $FLT_EPSILON))); if ($DEBUG) { print " last if (abs($phix) < $FLT_EPSILON && abs($phiy) < $FLT_EPSILON && abs($phiz) < $FLT_EPSILON && ($e == $e0 || abs($e+$e0) < $FLT_EPSILON || (abs(($e-$e0)/($e+$e0)) < $FLT_EPSILON))); \n"; } } #/*** Output data preparation ***/ for ($i=0; $i<$n; $i++) { $mol1->[$i+1][1] = $x1[$i]; $mol1->[$i+1][2] = $y1[$i]; $mol1->[$i+1][3] = $z1[$i]; } if ($DEBUG) { for ($i=0; $i<$n; $i++) { printf "%2.0f %10f %10f %10f\n", $m[$i], $x1[$i], $y1[$i], $z1[$i]; } printf "2*energy: %g\n", $e; } #/*** Finished ***/ # return sqrt($e); return $e; } sub proections { ############################################################################ ## Принимает ссылкy на молекулу и вторым параметром - плоскость, ## на которую делать проекцию ('xy','xz' или 'yz'). ## Возвращает сортированный по длинам проекций массив ## [проекция на плоскость xy вектора 0,0,0-атом , номер атома] ############################################################################ my $mol = shift; my $plane = shift; $plane ||= 'xy'; my @d; for (my $i=1; $i<@$mol; $i++) { my $d; if ($plane eq 'xy') { $d = sqrt($mol->[$i][1]**2 + $mol->[$i][2]**2); } elsif ($plane eq 'xz') { $d = sqrt($mol->[$i][1]**2 + $mol->[$i][3]**2); } elsif ($plane eq 'yz') { $d = sqrt($mol->[$i][2]**2 + $mol->[$i][3]**2); } push @d, [$d,$i]; } @d = sort {$a->[0] <=> $b->[0]} @d; return @d; } sub symmetr { ############################################################################ ## Принимает ссылку на молекулу. ## В свойства дописывает элементы хэша Symmetry => 'группа симметрии' ## Если определен 2-й параметр, то координаты будут симметризованы ## Требуется symmetry. ############################################################################ my ($mol,$symm_coord) = @_; my $N = @$mol - 1; open TINP, ">", "symm.TINP" or warn "Can't write temporary file: $!\n", return undef; print TINP " $N\n"; for (my $i=1; $i<=$N; $i++) { print TINP element($mol->[$i][0]), " $mol->[$i][1] $mol->[$i][2] $mol->[$i][3]\n"; } close TINP; my $prim = 0.1; my $final = 0.09; my @arg = ($symmetry_x, '-na', '-primary', $prim, '-final', $final, 'symm.TINP'); open OUT, '-|', @arg or do {warn "Can't run $symmetry_x: $!\n"; return undef}; while (<OUT>) { if (m/^It seems to be the (\w+) point group$/) { $mol->[0]{'Symmetry'} = $1; } if ($symm_coord && m/^Symmetrized coordinates of all atoms/) { <OUT>; my $nn = <OUT>; chomp $nn; for (my $i=1; $i<=$nn; $i++) { my ($nat,$x,$y,$z) = split ' ', <OUT>; print "$nat,$x,$y,$z\n"; $mol->[$i][0] = element($nat); $mol->[$i][1] = $x; $mol->[$i][2] = $y; $mol->[$i][3] = $z; } } } close OUT; $mol->[0]{'Symmetry'} ||= 'C1'; unlink 'symm.TINP'; return $mol->[0]{'Symmetry'}; } sub obminimize2xyz { my $mol = shift; my $opt_mol; if ($mol->[0]{File_name} =~ /\.sdf/) { my $args = "$obminimize_x -ff $method $MM_key -oxyz $mol->[0]{File_name} 1> _from_obmin_.xyz 2> _from_obmin_.out"; #print "$args\n"; system($args); } else { write_molden($mol, '_for_obmin_.xyz'); my $args = "$obminimize_x -ff $method $MM_key -oxyz _for_obmin_.xyz 1> _from_obmin_.xyz 2> _from_obmin_.out"; #print "$args\n"; system($args); #system($args) == 0 or warn "Err$? "; unlink('_for_obmin_.xyz'); } ($opt_mol) = read_molden('_from_obmin_.xyz'); return undef if @$opt_mol<2; open L, '<', "_from_obmin_.out" or die "Can't open obminimize output: $!\n"; while (<L>) { if (m/^\s*\d+\s+(-?\d+\.\d+)\s+-?\d+\.\d+$/) { #print "$1\n"; $opt_mol->[0]{'Energy'} = $1; } } close L; $opt_mol->[0]{'Point'} = $mol->[0]{'Point'} if $mol->[0]{'Point'}; unlink('_from_obmin_.xyz', '_from_obmin_.out'); return $opt_mol; } sub tminimize2xyz { my ($mol,$tinker_mm,$tinker_rms) = @_; return undef unless $mol; $tinker_mm ||= "$tinker_params/mmff"; $tinker_rms ||= 0.01; my @opt_mol = ($mol->[0]); return undef unless -x $babel_x; return undef unless -x $sdf2xyz2sdf_x; my $debug = 0; if ($mol->[0]{File_name} =~ /\.sdf$/) { open IN, '<', $mol->[0]{File_name} or die "Can't open $mol->[0]{File_name}: $!\n"; open OUT, '>', "_for_tmin_.sdf" or die "Can't write to _for_tmin_.sdf: $!\n"; while (<IN>) { print OUT $.==1 ? "_for_tmin_\n" : "$_"; } close OUT; close IN; } else { write_molden($mol, '_for_tmin.xyz'); open STDERR, '>>', "_for_tmin_.err" or die "Can't open _for_tmin_.err: $!\n"; my @args = ($babel_x, '--title','_for_tmin_', "_for_tmin.xyz", "_for_tmin_.sdf"); print "Run babel to covert xyz into sdf\n" if $debug; system(@args) == 0 or do {warn "System\n@args\nfailed: $?\n"; return undef}; } my $arg = "$sdf2xyz2sdf_x < _for_tmin_.sdf"; print "Run sdf2xyz2sdf to convert sdf into txyz\n" if $debug; system($arg) == 0 or do {warn "System\n$arg\nfailed: $?\n"; return undef}; rename '_for_tmin_.xyz', '_for_tmin_.txyz'; # Workaround for mmff-pibond bug if sdf2tinker open L, '<', "_for_tmin_.key" or die "Can't open _for_tmin_.key: $!\n"; open LL, '>', "_for_tmin_.key.tmp" or die "Can't write to _for_tmin_.key.tmp: $!\n"; while (<L>) { if (/mmff-pibond/i) { my @nn; while (/\s+(\d+)\s+(\d+)/g) { push @nn, $1<$2 ? [$1,$2] : [$2,$1]; } @nn = sort {$a->[0] <=> $b->[0] || $a->[1] <=> $b->[1]} @nn; print LL 'mmff-pibond'; foreach my $nn (@nn) { printf LL "%6d%6d", $nn->[0], $nn->[1]; } print LL "\n"; } else { print LL; } } close LL; close L; rename "_for_tmin_.key.tmp", "_for_tmin_.key"; my @args = ($tminimize_x, "_for_tmin_.txyz", $tinker_mm); push @args, 'A','A' if $tminimize_x =~ /newton$/; push @args, $tinker_rms; #system(@args) == 0 or do {warn "System\n@args\nfailed: $?\n"; return undef}; #print "@args\n"; open TINKER_OUT, '>', "$mol->[0]{File_name}.out" or die "Can't open $mol->[0]{File_name}.out: $!\n"; open TINKER, '-|', @args or do {warn "Can't run @args : $!\n"; return undef}; while (<TINKER>) { print TINKER_OUT; if (m/Final Function Value :\s*(\S+)/) { $opt_mol[0]{Energy} = $1; } } close TINKER; close TINKER_OUT; my $txyz = "_for_tmin_.xyz"; #my ($txyz) = map {$_->[0]} sort {$a->[1] <=> $b->[1]} map {[$_,-M]} glob("_for_tmin_.xyz*"); #print "$txyz\n"; open L, '<', $txyz or die "Can't open tminimize output $txyz: $!\n"; <L>; for (my $i=1; $i<@$mol; $i++) { my $line = <L>; my ($x,$y,$z) = (split ' ', $line)[2..4]; $opt_mol[$i][0] = $mol->[$i][0]; $opt_mol[$i][1] = $x; $opt_mol[$i][2] = $y; $opt_mol[$i][3] = $z; } close L; close STDERR; rename "_for_tmin_.sdf", "$mol->[0]{File_name}.sdf" if -f "_for_tmin_.sdf"; rename "_for_tmin_.err", "$mol->[0]{File_name}.err" if -f "_for_tmin_.err"; rename "_for_tmin_.key", "$mol->[0]{File_name}.key" if -f "_for_tmin_.key"; rename "_for_tmin_.xyz", "$mol->[0]{File_name}.1.txyz" if -f "_for_tmin_.xyz"; rename "_for_tmin_.txyz", "$mol->[0]{File_name}.0.txyz" if -f "_for_tmin_.txyz"; #print "$mol->[0]{File_name}\n"; unless ($debug) { foreach my $ext (qw(sdf err key 0.txyz 1.txyz out)) { unlink "$mol->[0]{File_name}.$ext" if -f "$mol->[0]{File_name}.$ext"; } } unlink "_for_tmin.xyz"; return \@opt_mol; } sub txyz2xyz { #my %mmff_atoms = ( #'CR'=>'C', 'HC'=>'H', 'OR'=>'O', 'NR'=>'N', 'S'=>'S', #'C=C'=>'C', 'HSI'=>'H', 'OC=O'=>'O', 'N=C'=>'N', 'S=C'=>'S', #'CSP2'=>'C', 'HOR'=>'H', 'OC=C'=>'O', 'N=N'=>'N', 'S=O'=>'S', #'C=O'=>'C', 'HO'=>'H', 'OC=N'=>'O', 'NC=O'=>'N', '>S=N'=>'S', #'C=N'=>'C', 'HOM'=>'H', 'OC=S'=>'O', 'NC=S'=>'N', 'SO2'=>'S', #'CGD'=>'C', 'HNR'=>'H', 'ONO2'=>'O', 'NN=C'=>'N', 'SO2N'=>'S', #'C=OR'=>'C', 'H3N'=>'H', 'ON=O'=>'O', 'NN=N'=>'N', 'SO3'=>'S', #'C=ON'=>'C', 'HPYL'=>'H', 'OSO3'=>'O', 'NR+'=>'N', 'SO4'=>'S', #'CONN'=>'C', 'HNOX'=>'H', 'OSO2'=>'O', 'NPYD'=>'N', '=SO2'=>'S', #'COO'=>'C', 'HNM'=>'H', 'OSO'=>'O', 'NPYL'=>'N', 'SNO'=>'S', #'COON'=>'C', 'HN'=>'H', 'OS=O'=>'O', 'NC=C'=>'N', 'STHI'=>'S', #'COOO'=>'C', 'HOCO'=>'H', '-OS'=>'O', 'NC=N'=>'N', 'S-P'=>'S', #'C=OS'=>'C', 'HOP'=>'H', 'OPO3'=>'O', 'NC=P'=>'N', 'S2CM'=>'S', #'C=S'=>'C', 'HN=N'=>'H', 'OPO2'=>'O', 'NC%C'=>'N', 'SM'=>'S', #'C=SN'=>'C', 'HN=C'=>'H', 'OPO'=>'O', 'NSP'=>'N', 'SSMO'=>'S', #'CSO2'=>'C', 'HNCO'=>'H', '-OP'=>'O', 'NSO2'=>'N', 'SO2M'=>'S', #'CS=O'=>'C', 'HNCS'=>'H', '-O-'=>'O', 'NSO3'=>'N', 'SSOM'=>'S', #'CSS'=>'C', 'HNCC'=>'H', 'O=C'=>'O', 'NPO2'=>'N', '=S=O'=>'S', #'C=P'=>'C', 'HNCN'=>'H', 'O=CN'=>'O', 'NPO3'=>'N', 'PO4'=>'P', #'CSP'=>'C', 'HNNC'=>'H', 'O=CR'=>'O', 'NC%N'=>'N', 'PO3'=>'P', #'=C='=>'C', 'HNNN'=>'H', 'O=CO'=>'O', 'NO2'=>'N', 'PO2'=>'P', #'CR4R'=>'C', 'HNSO'=>'H', 'O=N'=>'O', 'NO3'=>'N', 'PO'=>'P', #'CR3R'=>'C', 'HNPO'=>'H', 'O=S'=>'O', 'N=O'=>'N', 'PTET'=>'P', #'CE4R'=>'C', 'HNC%'=>'H', 'O=S='=>'O', 'NAZT'=>'N', 'P'=>'P', #'CB'=>'C', 'HSP2'=>'H', 'O2CM'=>'O', 'NSO'=>'N', '-P=C'=>'P', #'C%'=>'C', 'HOCC'=>'H', 'OXN'=>'O', '=N='=>'N', 'FE+2'=>'Fe', #'CGD+'=>'C', 'HOCN'=>'H', 'O2N'=>'O', 'N+=C'=>'N', 'FE+3'=>'Fe', #'CNN+'=>'C', 'HOH'=>'H', 'O2NO'=>'O', 'N+=N'=>'N', 'LI+'=>'Li', #'C5A'=>'C', 'HNR+'=>'H', 'O3N'=>'O', 'NCN+'=>'N', 'NA+'=>'Na', #'C5B'=>'C', 'HIM+'=>'H', 'O-S'=>'O', 'NGD+'=>'N', 'K+'=>'K', #'CO2M'=>'C', 'HPD+'=>'H', 'O2S'=>'O', 'NPD+'=>'N', 'ZINC'=>'Zn', #'CS2M'=>'C', 'HNN+'=>'H', 'O3S'=>'O', 'NR%'=>'N', 'ZN+2'=>'Zn', #'C5'=>'C', 'HNC+'=>'H', 'O4S'=>'O', 'NM'=>'N', 'CA+2'=>'Ca', #'CIM+'=>'C', 'HGD+'=>'H', 'OSMS'=>'O', 'N5A'=>'N', 'CU+1'=>'Cu', #'HN5+'=>'H', 'OP'=>'O', 'N5B'=>'N', 'CU+2'=>'Cu', #'HOS'=>'H', 'O2P'=>'O', 'N2OX'=>'N', 'MG+2'=>'Mg', #'HS'=>'H', 'O3P'=>'O', 'N3OX'=>'N', 'F'=>'F', #'HS=N'=>'H', 'O4P'=>'O', 'NPOX'=>'N', 'CL'=>'Cl', #'HP'=>'H', 'O4CL'=>'O', 'N5M'=>'N', 'BR'=>'Br', #'HO+'=>'H', 'OM'=>'O', 'N5'=>'N', 'I'=>'I', #'HO=+'=>'H', 'OM2'=>'O', 'NIM+'=>'N', 'F-'=>'F', #'OH2'=>'O', 'N5A+'=>'N', 'CL-'=>'Cl', #'OFUR'=>'O', 'N5B+'=>'N', 'BR-'=>'Br', #'O+'=>'O', 'N5+'=>'N', 'CLO4'=>'Cl', #'O=+'=>'O', 'N5AX'=>'N', 'SI'=>'Si', #'N5BX'=>'N', #'N5OX'=>'N', #'NPY'=>'N'); my %mmff_atoms = map {$_=>1} qw(C H O N S P F CL BR I LI NA K ZN CA CU MG SI); my ($file,$energy) = @_; my $mol; open F, '<', $file or do {warn "Can't open $file: $!\n"; return}; $mol->[0]{'Energy'} = $energy; <F>; while (<F>) { my ($at,$x,$y,$z) = (split)[1..4]; my $atom = $at; $at =~ s/^[^A-Z]+//; chop $at until exists $mmff_atoms{$at}; do {warn "Undefined MMFF atom $atom\n"; return} unless $at; push @$mol, [ucfirst(lc $at),$x,$y,$z]; } close F; return $mol; } sub sdf2mol { my $file = shift; my @mols; my $num = qr/^-?\d+(?:\.\d+)?/; open L, '<', $file or die "Can't open $file: $!\n"; while (<L>) { my $mol; my $charge = 0; my $id = $_; chomp $id; my $vendor = <L>; chomp $vendor; my $comment = <L>; chomp $comment; my $count_str = <L>; my ($N) = $count_str =~ /^\s*(\d+)/; die "Not SDF $file? ($count_str)\n" unless $N; for (my $i=1; $i<=$N; $i++) { my $str = <L>; #print $str; my ($x,$y,$z,$at) = map {s/\s+//g; $_} unpack("a10 a10 a10 a4", $str); #print "$x,$y,$z,$at\n"; if ($x !~ /$num/ || $y !~ /$num/ || $z !~ /$num/ || $at !~ /^\w?\w$/) { die "Not SDF $file? ($x,$y,$z,$at)\n"; } $mol->[$i] = [$at,$x,$y,$z]; } while (<L>) { if (s/^M\s+CHG\s+\d+\s+//i) { my %hhh = split ' ', $_; foreach my $ch (values %hhh) { $charge += $ch; } $mol->[0]{Charge} = $charge if $charge != 0; } if (/^\$\$\$\$/) { push @mols, $mol; last; } } } close L; #dd @mols; exit; return @mols; } sub rms_ellips { my @e1 = split /_/, $_[0]; my @e2 = split /_/, $_[1]; return sqrt(($e1[0]-$e2[0])**2+($e1[1]-$e2[1])**2+($e1[2]-$e2[2])**2); } ####################### # Иерархия функций ####################### #&main # &read_molden # &tenzor # &jacobi # &sum # %massa # &get_rms # &xy_proections # &small_RMS # &tenzor_sort # &znak # &kinen # &sum # %massa # &delete_same # &print_write_molden # &write_molde