Расчёт энергии сольватации димера из двух мономеров полиэтиленгликоля методом FEP

  1. Разобьем 1,2-диметоксиэтан на два мономера. SMILES исходного вещества: СОССОС. Мономер: СОС.
  2. Расчитаем частичные заряды на атомах этого остатка используя как заглушку метильную группу.
    SMILES: СОСC.
    obgen 1.smi > 1.mol
1.pdb
export PATH=${PATH}:/home/preps/golovin/progs/bin
Ante-RED.pl 1.pdb
    В полученный файл 1-out.p2n переименован в Mol_red1.p2n, где добавлено следующее:
    REMARK INTRA-MCC 0.0 | 9 10 11 12 | R
    Запущен скрипт: 
    RED-vIII.4.pl Mol_red1.p2n
  3. Построим запись для остатка (rtp) в силовом поле amber99sb.
    cp -r /usr/share/gromacs/top/amber99sb.ff/  .
    (копирование директории силового поля в рабочую директорию)
    Запись для глицина в файле amber99sb.ff/aminoacids.rtp:
    [ GLY ]
 [ atoms ]
     N    N           -0.41570     1
     H    H            0.27190     2
    CA    CT          -0.02520     3
   HA1    H1           0.06980     4
   HA2    H1           0.06980     5
     C    C            0.59730     6
     O    O           -0.56790     7
 [ bonds ]
     N     H
     N    CA
    CA   HA1
    CA   HA2
    CA     C
     C     O
    -C     N
 [ impropers ]
    -C    CA     N     H
    CA    +N     C     O
    Изменим файл aminoacids.rtp, добавив в конец запись нашего мономера (назовём его EGL):
    [ EGL ]
 [ atoms ]
    CA    CT           0.02750     1
   HA1    H1           0.05010     2
   HA2    H1           0.05010     3
   HA3    H1           0.05010     4
     O     O          -0.36850     5
    CB    CT           0.13770     6
   HB1    H1           0.02000     7
   HB2    H1           0.02000     8   
   
 [ bonds ]
    CA   HA1
    CA   HA2
    CA   HA3
    CA     O
     O    CB 
    CB   HB1
    CB   HB2
    CA   +CB
    Заряды были взяты из Mol_m1_o1-sm.mol2.
  4. Подготовка системы к МД и запуск счёта с шагом lambda 0.1.
    obgen 2.smi > 2.mol
babel -imol 2.mol -opdb 2.pdb
Ante-RED.pl 2.pdb
    В полученном PDB-файле был переименован остаток (LIG -> EGL), атомы названы также, как и в aminoacids.rtp. Поскольку получен PDB-файл для димера, то остаткам (мономерам) даны соответствующие номера. Подготовка системы к МД в пакете Gromacs: 
    cd lynx/md                                                         ##создана папка md* в корневой директории
export X3DNA=/home/preps/golovin/progs/X3DNA
export PATH=/home/preps/golovin/progs/X3DNA/bin:${PATH}
pdb2gmx -f 2-out.pdb -o dimer -p dimer -ff amber99sb -water tip3p  ##построим файл топологии системы в силовом поле amber99sb 
                                                                 и файл с координатами в формате Gromacs 
                                                                (структура димера находится в файле 2-out.pdb) 
editconf -f dimer.gro -o dimer_ec -d 1                             ##при построении ячейки сделан отсуп на 1 нм
grompp -f em -c dimer_ec -p dimer -o dimer_em -maxwarn 1           ##оптимизация геометрии системы для удаления "плохих" контактов в молекуле
mdrun -deffnm dimer_em -v                                          
genbox -cp dimer_em -p dimer -cs -o dimer_s                        ##добавление в ячейку молекулы воды
grompp -f em -p dimer -c dimer_s -o dimer_s                        ##нейтрализация заряда системы (genion не запускаем)
grompp -f pr -c dimer_s -p dimer -o dimer_pr -maxwarn 1            ##"утряска" воды
mdrun -deffnm dimer_pr -v
    Kороткое md на 100 пс (это время уже прописано в файле md.mdp) 
    grompp -f md -c dimer_pr -p dimer -o dimer_md -maxwarn 1
mdrun -deffnm dimer_md -v
    *в папке md: директория amber99sb.ff, 2-out.pdb, em.mdp, pr.mdp, md _lam.mdp, md_out.mdp.
    Cкрипт для запуска 10 траекторий с разным значением lambda. 
    g_energy -f lamda.edr -o lamda01.xvg
    Gnuplot: 



    trjconv -f lamda.trr -s lamda.tpr -o lamda01.pdb


  5. Определение значение энергии и оценка потенциально не точно посчитанных этапов. 
     g_bar -f */lamda.xvg -o -oi -oh                                   ##oпределение энергии перехода из димера в ничто
    Построение зависимости из файлов bar.xvg barint.xvg:
    Энергия гидратации не совпадает с экспериментальной 4.8 Kcal/mol. Для этого надо изменить параметры следующим образом:
© Anastasia Maslova, 2012