На главную

Вычисление точечных зарядов и VdW параметров для молекулярной механики

  1. .......

    2. при помощи babel получили pdb файл этана из результатов оптимизации из предыдущего практикума et.log.
    добавляем путь к скриптам: export PATH=${PATH}:/home/preps/golovin/progs/bin
    C помощью скрипта Ante_RED.pl подготовим pdb файл (Ante_RED.pl et.pdb). В файле et-out.p2n находим информацию о заряде и мультиплетности (0 и 1, соответственно. Все верно => переходим далее). Файл et-out.p2n переименовываем в Mol_red1.p2n. Запускаем RED (RED-vIII.4.pl) В результате работы скрипта получаем Mol_m1-o1.mol2 файл, содержащий информацию о координатах и зарядах атомов (и помимо этого файла много других).
    Создаем файл описание молекулы этана в формате пакета программ GROMACS et.top соуществляем его "размножение" 
    
for i in {1..7};do
    ep=$( echo "scale=5; 1/$i/$i/$i"  | bc -l )
    sed "s/тут надо что-то написать/$ep/" et.top > v_${i}.top
done

    Мы создали 7 файлов топологии. Теперь надо провести для каждой системы молекулярную динамику с каждым файлом топологии. Добавляем в скрипт строчки для расчета. 
     grompp_d -f md -c box_big -p v_${i}.top -o vb_${i} -maxwarn 1 && mdrun_d -deffnm  vb_${i} -v 
 grompp_d -f md -c box_38 -p v_${i}.top -o v_${i} -maxwarn 1 && mdrun_d -deffnm  v_${i} -v

    Проведем расчет и комментируем строчки со счетом в скрипте. Можно конвертировать полученные trr файлы в pdb и посмотреть на них в PyMol. Для конвертирования воспользуемся следующей командой:
    trjconv_d -f v_1 -s v_1 -o v_1.pdb
    получаем следующие файлы:
    v_1.pdb
    v_3.pdb
    vb_1.pdb
    vb_3.pdb
    v_1 vs v_3 (жидкий этан).
    Исходно все молекулы сконцентрированы в одном месте. при запуске анимации молекулы начинают распределяться в пространстве (для всех файлов, кроме v_1). Некоторые молекулы перемещаются группами. в случае v_1 видим, что небольшая часть молекул покинула свое исходное положение. все же остальные остались в исходной точке. vb_1 vs vb_3 (газообразный этан)
    Изначально молекулы не сконцентрированы в определенной точке пространства. после запуска анимации молекулы начинают двигаться (при этом есть доля молекул, которые практически не изменяют своего положения). между собой случаии не отличаются. молекулы движутся одинаково. И лишь одна молекула практически не меняет сввоего положения.
    Теперь нам надо посчитать сами значения энергий, для этого воспользуемся утилитой g_energy. Эта утилита может работать в интерактивном режиме, но это не удобно в скрипте поэтому используем перенаправление потока. Символ '\n' означает перенос строки:
    echo -e "LJ-(SR)\nCoulomb-(SR)\n0" | g_energy -f -b 10 vb_${i} -o eb_${i} > vb_${i}.txt echo -e "LJ-(SR)\nCoulomb-(SR)\n0" | g_energy -f -b 10 v_${i} -o e_${i} > v_${i}.txt
    В результате получаем следующие зависимости энергии 

    		Жидкий этан			Газообразный	
epsilon		LJ (kJ/mol)	Coulomb(kJ/mol) LJ(kJ/mol)	Coulomb (kJ/mol)
1.00000  	-253.389	0.00253608	-0.000226538	3.89948e-05  
0.12500  	-31.9072 	0.000452278	-0.000157918	3.89893e-05  
0.03703	 	-5.14133 	0.000149597	-0.000141251	3.89911e-05  
0.01562	 	-1.57194 	8.86437e-05	-0.000134331	3.89893e-05  
0.00800 	-3.99826 	2.56685e-05	-0.000130696	3.89908e-05   
0.00462  	-2.60282 	0.000114806	-0.000128504	3.89909e-05   
0.00291  	-3.44165 	4.65368e-05	-0.000127074	3.89911e-05
    В жидкой фазе наибольший вклад вносят VdV взаимодействия и Кулоновские взаимодействия пренебрижимо малы. в случае газообразной фазы - энергия VdV взаимодействий отличается от кулоновских лишь на порядок. Энергии VdV взаимодействий значительны при низких температурах (жидкое состояние), а в газообразно фазе их значения существенно снижаются.
    Вклад VdV взаимодействий в газообразной фазе и Кулоновских взаимодействий, как в газообразной, так и в жидкой фразе, в энергию испарения этана незначителен, следовательно пренебрежем ими.
    Значение epsilon может находится в интервале от 0.03703 до 0.12500
    Точное значение epsilon
    Для вычисления точного значения epsilon модифицируем скрипт (а именно первые 2 строчки):
    for i in {1..10};do
    ep=$( echo "scale=5; 0.03703+($i-1)*0.005" | bc -l )
    с таким шагом не получил удовлетворительного результата. пробуем еще поменять шаг.
    уменьшим шаг и повторим операцию (
    for i in {1..5};do
    ep=$( echo "scale=5; 0.03703+($i-1)*0.0005" | bc -l )     )
    после нескольких итреаций изменения скрипта удалось вычислить наиболее оптимальное значение epsilon, при котором энергия испарения соответствует экспериментальным данным. это достигается при epsilon(для водорода)=0.04454. Энергия испарения при этом значении эпсилон равно 5.45kJ/mol
    ©Анисенко Андрей