Transcript Энергия сольватации

Энергия сольватации
соответствие между гидродинамическим и
электростатическим описанием
•
•
•
•
•
+4πQ -положительный точечный заряд,
умноженный на 4π – какая масса жидкости
вытекает из точки за единицу времени
- 4πQ -отрицательный точечный заряд ,
умноженный на 4π – какая масса жидкости втекает
в точки за единицу времени
ε - диэлектрическая проницаемость – плотность
жидкости (При использовании диэлектрической
проницаемости объемный и поверхностный заряд
диэлектрика) Для ваакума плотность жидкости
равна единице
E - напряженность электрического поля – скорость
жидкости
∫EdS - поток поля через поверхность – какая масса
жидкости протекает через поверхность за единицу
времени
Физический смысл теорема
Остроградского –Гаусса
• Поток поля через замкнутую
поверхность равен сумме зарядов
внутри поверхности, умноженной на 4π
- имеет теперь тривиальный смысл
закона сохранения массы жидкости!
Сколько жидкости вытекает из
источников в сумме, столько же
вытекает за поверхность.
• ∫SEdS=4π∑jQj
Физический смысл точечного
заряда – закон Кулона
• ε E S= 4πQ  ε E 4πR2= 4πQ
Отсюда закон Кулона:
• E=Q/ (ε R2 )
Нормальная составляющая
поля на границе диэлектрика
• ε1E1n dS= ε2E2n dS
• Отсюда соотношение для нормальных
компонент поля
• ε1E1n= ε2E2n
Поверхностный заряд (без замены
поверхностного заряда на ε):
Поле плоской заряженной
поверхности в вакууме
• 2 E S=4πσ S
• E=2πσ
Потенциал φ

 
Ed l  0
E = -grad φ
 
 
( rot E  n )  Ed l / dS  0

Условие для тангенциальной
составляющей на границе диэлектрика
• Еτ1∆l – Еτ2∆l =0
• Еτ1=Еτ2
Уравнение Лапласа

 
EdS  Q  dV
4 

 
EdS / dV 

 div( εE)   div( ε grad  )
Разбиение поверхности на внутренние
поверхностные элементы
Введем обозначения
• qк - заряд к-ого элемента.
• Sк – площадь поверхностного элемента
• rj – координата центра поверхностного
элемента
• nj – нормаль в центре поверхностного
элемента
• Qi – заряды внутри молекулы
• Rj – координата зарядa j внутри
молекулы
Метод РСМ- поле диэлектрика
описывается поверхностным зарядом
полости
• Поток поля поверхного элемента с
номером k≠m через поверхностный
элемент с номером m
• (Е1n)km Sm = (qk (rkm∙nm) /(rkm3) Sm
rkm= rm- rk
Суммарный поток поля всех точечных
зарядов внутри полости через
поверхностный элемент m


( E1n )Q Sm  



i

  
Qi ( rm  R i ) n m 
S

m

3


in rm  R i



Поток собственного поля
элемента m
(E1n)m,self Sm = -2π (qm/ Sm) Sm=-2π qm
(E1n)m,self Sm = -2πqm[1-(Sm/{4πRm2})1/2 ]
∑ k≠m ((Е1n)mk Sk) + (E1n)m,self Sm =0
(E1n)m,self Sm = - ∑ k≠m ((Е1n)mk Sk) =
=- qm ∑ k≠m (((rk- rm) ∙nk) / | rk- rm |3)Sk
Полный поток поля
• qm = (1- εin/ εout) E1n Sm/(4π) 
• E1n Sm=-4πqm/ (1- εin/ εout)
• С другой стороны он равен сумме всех
его трех составляющих
• E1n Sm = (∑k≠m (Е1n)km) Sm  ( E1n ) Q Sm 
+ (E1n)m,self Sm
Вывод уравнения РСМ
уравнение РСМ в матричной
форме
A(PCM) q = -B(PCM) Q
Вывод уравнения COSMO
• Потенцилы всех поверхностных элементов в точке
поверхностного элемента m, потенциал всех
зарядов внутри полости на поверхностном
элементе m и собственный потенциал
поверхностного элемента m в сумме равны нулю.
φmm - Собственный потенциал поверхностного элемента
в матричной форме
Собственный потенциал
поверхностного элемента
mm  qm  3.8 | Sm
1

| 2
Из условия 1)
 ni Si  q j
Norm
j i
r  r   q
i
ri  rj
j
3
Norm
i
 r  r 

j
i

njS j   0

3
 j i r  r

j
i


Условия нормализации
СОSMO.
Поверхность с укрупненными
поверхностными элементами
Нулевое приближение. Заряд в идеально
проводящей сферической полости. Точное
решение методом зеркального заряда.

1
1
RQ 2
1 Q2
1
1 Q2
G  Q  ( rs )  


2
2
2
2 R
2 R 
2 a


r
s
1 

rs 
 rs 
2
 rs

 R 




a  R (1 
G

1
2

1
2

i
j

i
1
2

1
Qi j ( ri )  
2
QiQ j
j

i
j
QiQ j

i
j
 
( rjri )  R  2rjri R cos(r1, r2 )
1
R ij
4
R2
)
R
1
QiQ j 

| rj |  R 2 
ri   2 rj
rj
R
2
rs2
2

1
2

i
j
QiQ j
1
 
( rj  ri ) 2  a ia j
Уточнение первое.Учет диэлектрической
проницаемости
2
Q
1
1
G   1  
2
 a
1 1
G   1  
2 
n
n
 R
Qi Q j
i 1 j1
1
ij
1
 n n Qi Q j
1 1
E  

2  in 1  1 i 1 j 1 Rij
2
Учет ионной силы раствора
 Rij

1 1  e
G  
1

2  in


1 1
G  
2 in
n
 e0.73R ij
1 


j1 
n

i 1
Q 2

 a

 QiQ j

 R ij

Уточнение 3. Учет
несферичности полости
поверхностный заряд
1
1
(r )   (r ) 
 Qj  
 in i
r  rj
SGB
SES и SAS
SES (Solvent Excluded Surface) - поверхность исключённого из
растворителя объёма.. Объем, занимаемый растворителем лежит
вне объема, ограниченного этой поверхностью. Сам субстрат
полностью лежит внутри этого объема.
SAS (Solvent Accessible Surface) - поверхность доступная
растворителю образуется центрами молекул растворителя,
касающихся молекулы субстрата.
Электостатическая часть должна описываться SES, поскольку
SAS плохо описывает многие уччастки молекулы, существенные
для взаимодействия [6].
Поверхность SES молекулы можно описывать [1b] (1) гладко,
заменяя ее простыми формами типа сферы, элипсоида или
цилиндра, (2) детально (а) покрытием из ван-дер-ваальсовских
сфер вокруг атомов (b) покрытием из сфер вокруг химических
групп атомов (с) как в предыдущих двух методах, но заполняя
остающееся пустое пространство внутри SES фиктивными
сферами (GEPOL) [57-58] (d) соединяя сферы, описанные в (а)
или (b) участкам вогнутых поверхностей [59-60].
Преобразование поверхности типа SES в поверхность типа SAS
Неэлектростатическая часть
энергии в модели Абагяна
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
Пусть имеется поверхность после первичной обкатки.
Пусть имеется треугольник на поверхности SES (r1,r2,r3)
Образ этого треугольника на SАS (r1n,r2n,r3n)
r1n=r1+n1*p_rol
(1)
r2n=r2+n2*p_rol
(2)
r3n=r3+n3*p_rol
(3)
где p_rol – радиус первичной обкатки,
n1, n2, n3-нормали в соответствущий точках.
Тогда
1) поверхностный элемент на первичной сфере, опирающейся на три атома, отображается в точку. Его
площадь на SАS нулевая
2) поверхностный элемент на торе, опирающимся на два атома, отображается в линию. Его площадь на
SАS нулевая
3) поверхностный элемент на атоме отображается в сферический тругольник. Его площадь на SАS
пропорцианальна его площади [(ratom+ p_rol)/ ratom]2 Градиент площади по смещению атомов равен нулю
4) граничный поверхностный элемент между тором и первичной сферой, опирающейся на три атома –
нулевая площадь и градиент
5) граничный треугольный поверхностный элемент между тором (и/или первичной сфере, опирающейся на
три атома) и атомом – площадь считается по формуле (187) для трех точек (r1n,r2n,r3n)
ESAS=∑j σj sj + b (1)
j=1,N по всем (ненулевым на SAS) поверхностным элементам
Для воды:
σj =σ =0.00378 - в ккал/(моль А2)
b=0.698
- в ккал/моль
Вторичная обкатка поверхности
Применение метода вторичной обкатки для
сложной геометрической конфигурации
расположения атомов
Точки на атомах
Ван-дер-Ваальсовская поверхность
Тройная точка
Входные данные: координаты и Ван-дер-ваальсовые радиусы атомов,
радиус сферы обкатки.
Цикл перебора
троек атомов c
индексами i,j,k;
i≠j, i≠k, j≠k
Вычисление координат центра сферы обкатки, касса-ющейся текущей тройки
атомов по формулам (340)-(349)
Проверка пересечения сферы обкатки с центром с каким-либо атомом молекулы, не в
ходящим в текущую тройку опорных атомов.
Запомнить координаты и
индексы текущих трёх опорных
атомов.
Да
Да
Обработаны все
тройки атомов?
Есть
пересечение?
Нет
Нет
Перейти к обработке
следующей тройки атомов.
Выход из цикла.
Выход из процедуры с выдачей массива координат центров вогнутых сферических элементов между тройками
атомов.
Первичная обкатка двух
атомов
Входные данные: координаты и радиусы атомов, координаты положений
центра сферы обкатки при контакте с тремя атомами.
Цикл перебора
пар атомов.
Вычисление координат центра и радиуса круговой траектории обкатки вокруг
текущей пары атомов. Задание вектора z локального базиса. (350)- (353)
Поиск атомов нарушающих свободную обкатку сферы вокруг данной пары атомов.
Определение
дополнительных
базисных векторов
согласно выражению
(353).
No
Найден хотя
бы один атом.
Да
Определение
дополнительных базисных
векторов из условий (354).
Определение координат центра сферы обкатки при касании трёх атомов, два из которых составляют
текущую пару, а третий – один их нарушающих свободную обкатку. Обработка всех конфликтных с
текущей круговой траекторией атомов.
Задание массива углов согласно (355) – для всех найденных положений центра сферы обкатки при
контакте с тремя атомами, анализ и определение начальных’’ и конечных ’’ углов для каждой
дуги соотнесённой к текущей паре атомов.
Запомнить вычисленные параметры в соответствующих ячейках массивов структур данных.
Да
Обработаны все
пары атомов?
Выход из цикла.
Нет
Переход к обработке
следующей пары
атомов.
Выход из процедуры с выдачей массивов структур данных описывающих параметры тороидальных
фрагментов.
Рисунок . Блок схема алгоритма определяющего
массив параметров тороидальных фрагментов.
Точки на поверхности вторичных
сфер устойчивых положений
вторичная обкатка для точек
находящихся на сфере шар-зонда при
его опоре на три первичные сферы


 ( r  r1 ) 2  a 2

 2
 
(
r

r
 b2

2)
 ( r  r ) 2  c 2
3


 

r  pc  zh
 
 

(r  r )  (r2  r1)
z  3 1
(r2  r1)  (r3  r1)


 ( pc  r1 ) 2  a 2  h 2

 2
 
2
2
( pc  r2 )  b  h

 2
( p

r
)  c2  h2
c
3


( r3  r1)  ( r2  r1 )  ( r2  r3 )  (r1  r3 )  (r1  r2 )  ( r3  r2 )
  





(
r r r )
h 2  (a 2  b2  c2 ) / 3  ( r12  r22  r32 ) / 3  pc2  2pc  1 2 3
3


 
(
r

r
 2
1) 
 (r  r ) 
1
 3



pc  (a 2  b 2  r12  r2 2 ) / 2



pc  (a 2  c 2  r12  r3 2 ) / 2
(1)


  
  
pc  r1  1[ z  ( r2  r1 )]  1[z  ( r3  r1 )]
 

a 2  c 2  ( r3  r1 ) 2
1 



 
2 ( r2  r1 )  ( r3  r1 )




a 2  b2  ( r2  r1 ) 2

1 


 


2
(
r

r
)

(
r3  r1 )
2
1





2
2

b  c  ( r3  r1 ) 2  ( r2  r1 ) 2
 1  1  1 


 

2
(
r

r
)

(
r3  r1 )

2
1


 



pc  r1  [z  ( 1r1 1 r2  1r3 )]
 2 2   2   2






3
(
b

c
)

(
r

r
)

(
r

r
)


3
1
2
1



[
(
r

r
)

(
r

r
)

r
]

3
1
2
1
1
    2




(
r

r
)

(
r

r
)
3
1
2
1


  
 
 


 ( r1  r2  r3 ) 1  3(c2  a 2 )  ( r1  r2 )2  ( r3  r2 )2      
pc 
 
[( r1  r2 )  ( r3  r2 ) r2 ]  




2


3
6

( r1  r2 )  ( r3  r2 )






 3(a 2  b2 )  ( r2  r3 )2  ( r1  r3 )2      
[( r2  r3 )  ( r1  r3 ) r3 ] 
 
    2


( r2  r3 )  ( r1  r3 )


Вторичная обкатка двух первичных
сфер.
Торы вторичной обкатки.
Определение положения затравочного треугольника
Положение затравочного треугольника определяет
начало построения сетки триангуляции. С точки
зрения алгоритма построения это положение может
быть выбрано произвольно. Задаётся некоторая
удалённая от молекулы точка в пространстве, для
которой определяется соответствующая проекция на
поверхность. Затем вокруг найденной точки проекции
выстраивается равносторонний треугольник с
размером радиуса описанной окружности примерно
равной величине заданной пользователем. Данные о
первом треугольнике помещаются в массивы
структур данных описывающих сетку триангуляции
Процедура обхода рёбер текущей замкнутой ломаной.
Цикл перебора текущих ребер и выбор подходящего сценария
(список всех возможных сценариев приведен ниже описания алгоритма):
1. Выбираем текущее ребро следующим образом:
Предположим что ранее, после добавления нового треугольника текущая замкнутая ломаная разделилась на
две новые замкнутые ломаны с одной общей вершиной. Тогда в качестве текущего ребра выбираем одно из
ребер, имеющих одним из концов общую точку этих двух образовавшихся ломаных. Целью такого выбора
является пространственное разделение двух ломаных. Действительно, добавление новых треугольников
приводит к исчезновению их общей точки. После образования двух не связанных ломаных прекращаем
выполняемый цикл обработки текущей ломаной и идем на начальный пункт перебора всех ломаных. Если
разделения ломанной на две не произошло, то берем следующее ребро после текущего в текущем множестве
ребер.
2.Пусть текущее множество, включающее текущее ребро, состоит из трех ребер. Тогда используем сценарий 10
(и идем на пункт 1).
3. Пусть текущее множество, включающее текущее ребро, состоит из четырех ребер.
Пусть у образованного ими четырехугольника имеется «особая» вершина, обладающая двумя свойствами. Вопервых, из нее выходит только два ребра. Во-вторых, на нее опирается только один треугольник, образованный
этими двумя ребрами. Если двухгранный угол между этой треугольной гранью с «особой» вершиной и гранью с
противоположной ей вершиной четырехугольника меньше π/6, то используем сценарий 9 (и идем на пункт 1), а
если больше – то сценарий 8 (и идем на пункт 1). Если «особой» точки нет, то используем сценарий 7 (и идем на
пункт 1).
4.Один из углов текущего ребра с одним из двух смежных ребер из текущего множества не маленький (α1> π/9), а
с другим маленький (α2< π/9):
Если разные концы текущего ребра и смежного ребра, образующие малый угол не связанны еще одним путем из
двух ребер, то используем сценарий 6 (и идем на пункт 1). Иначе используем сценарий 5 (и идем на пункт 1).
5.Оба угла α1 и α2 – маленькие (α1< π/9 и α2< π/9), то используем сценарий 1 (и идем на пункт 1).
6.Строим «новый» треугольник и «новую» точку по сценарию 2 , но пока не подтверждаем их построение в
качестве узла и треугольника сетки.
7.Пусть углы β1 и β2 - углы между ребрами, смежными к текущему ребру, и соответствующими смежными к ним
сторонами «нового» треугольника
8.Пусть хотя бы один из углов β1 или α1 мал, а углы β2 и α2 -велики. ([(β1<π/6 или α1<2π/9) и (β2>π/6 и α2>2π/9)]).
Или, наоборот, хотя бы один из углов β2 или α2 мал, а углы β1 и α1 -велики [(β2<π/6 или α2<2π/9) и (β1>π/6 и
α1>2π/9)]). Тогда применяем сценарий 4 для текущего ребра и смежного ребра, образующего малый угол (и идем
на пункт 1).
Продолжение алгоритма
12.
Пусть углы β1 и β2 малы (β1<π/6 и β2< π/6). Тогда применяем сценарий 4 для текущего ребра и смежного ребра, образующего
меньший угол (и идем на пункт 1).
. Пусть все определенные выше углы β1, α1, β2 и α2 велики (β1>π/6, α1>2π/9 , β2>π/6 и α2>2π/9). Пусть nmid - нормали к
поверхности в точке проекции середины текущего ребра, nnew – нормаль к поверхности в «новой» точке «нового» треугольника.
Пусть γ- угол между nnew и nmid. Если угол γ велик (γ>π/2) , то применяем сценарий 1 (и идем на пункт 1).
Пусть хотя бы один из определенных выше углов β1, α1, β2 и α2 мал (β1<π/6, или α1 <2π/9 , или β2<π/6, или α2<2π/9). Тогда
выполняем сценарий 1 (и идем на пункт 1).
Находим «первую особую» точку, близкую к двум вершинам текущего ребра по формуле:
,
r( 1 ) 
9.
10.
11.
r1  r2  2rnew
4
(6)
rnew – радиус-вектор «новой» точки.
r1, r2 - радиус-вектора двух вершин текущего ребра. Если для всех точек (в кубической области, содержащей «первую особую» точку)
Формируем множество всех узлов сетки, которые обладают следующими свойствами:
а. Лежат в кубической области, содержащей «первую особую» точку
в. Не совпадают с вершинами текущего ребра и вершинами двух смежных с ним граничных ребер
Находим «вторую особую» точку. Эта точка из найденного множества, имеющая минимальное расстояние ||rmin|| до «первой особой»
точки.
13. Если расстояние ||rmin||>Rch, то подтверждаем построение «нового» треугольника по сценарию 2 (и идем на пункт 1). Если
α1<π/2 или α2<π/2 , то выполняем сценарий 1 (и идем на пункт 1).
14. Если «вторая особая» точка не лежит на текущей границе, то выполняем сценарий 1 (и идем на пункт 1).
15. Если ||rmin|| меньше половины текущего ребра, то выполняем сценарий 1 (и идем на пункт 1).
16. Пусть nm – нормаль во «второй особой» точке. Если угол между nnew и nm больше π/2 , то выполняем сценарий 1 (и идем на
пункт 1).
17. Выполняем сценарий 3. Образуются две замкнутые ломаные граничных ребер с общей точкой. Идем на пункт 1.
Основные сценарии добавления нового треугольника или обработки текущих граничных рёбер:

R g1 R1g 2 R2g1 R2g 2 
L  m in Lmin , 1
,
,
,
(7)
2
2
2
2 


R g1 R g 2 
L  m in L , N , N  ( 8 )
2
2
 2
r
N a ra  N b rb 
N a  N b  ( 9 ),( 10 ),( 11 )
1)
2)
Ничего не делать. Когда все другие сценарии не сработали.
Строится «новый» треугольник, опирающийся на текущее ребро. (Рис. 2) Делается это следующим образом. Проецируем центр текущего ребра на поверхность.
Из этой точки строим вектор, перпендикулярный нормали в этой точке и вектору текущего ребра и длинной равный Rch =1.5L. Проецируем полученную точку на
поверхность. Это «новая» точка. Строим «новый» треугольник из текущего ребра и полученной «новой» точки.
L – это адаптивный шаг сетки. Адаптивный шаг сетки - это радиус окружности, описанной вокруг равностороннего треугольника с высотой Rch =1.5L. Этот
треугольник определяет максимальный размер триангуляции в данном месте.
Адаптация шага сетки идет по следующему алгоритму. Если к кубу, определяемого серединой текущего ребра, не относится ни один «центр адаптации», то этот
шаг сетки определяется заранее заданным максимальным размером L =Lmax. Если такие центры есть, то считается расстояние от каждого из этих «центров
адаптации», относящихся к кубу, до середины текущего ребра ||rj|| , (j=1,…,Na – номер центра адаптации). Находим среди этих расстояний те, которые
меньше критического для соответствующего «центра адаптации» ||rj|| < Rаj. Каждому такому «центру адаптации» j соответствует свой шаг сетки L=Lаj.
Выбираем среди них минимальный шаг L=Lmin.
Каждой двух вершин текущего ребра (впрочем, как и любой точке поверхности) соответствует два главных радиуса кривизны Rg1, Rg2. Для тора это радиусы двух
образующих тора в этой точке (один из них всюду одинаков, другой увеличивается от центра к краям тора). Для сферического сегмента оба главных радиуса
равны его радиусу. Корректируем шаг сетки, чтобы он был не больше половины этих радиусов для обеих вершин:
.
(7)
Далее строим «новую» точку с таким шагом сетки.
Пусть «новая» точка лежит на сегменте (тороидальном или сферическом) SN. Вершины текущего ребра лежат на сегментах S1, S2. Проверяем, что SN либо
совпадает с S1 или S2, либо является соседним для обоих из них. Для «новой» точки находим ее два главных радиуса кривизны Rg1N, Rg2N. Проверяем, что
L< Rg1N/2, L< Rg2N/2. Если хотя бы одно из этих условий не выполняется, то шаг сетки уменьшается следующим образом:
.
(8)
Далее строим «новую» точку с таким шагом сетки и снова проверяем описанные выше условия. Этот процесс продолжается до тех пор, пока эти условия не
выполнятся.
3)
«Новая» точка»-а, построенная по методу 2 «сливается» в одну точку со «второй особой» точкой (узлом сетки) - b. «Слияние» идет по формуле:
,
(9)
Nа =2 - число ребер, выходящих из точки а; Nb - число ребер, выходящих из точки b;
rа – радиус-вектор точки а; rb - радиус-вектор точки b; r - радиус-вектор образующейся в результате слияния точки;
4) Строится новый треугольник. Он образуется, во-первых, текущим ребром. Во-вторых, граничным ребром, смежным к текущему ребру и имеющим с ним малый угол.
И, в-третьих, одним новым граничным ребром, построенным напротив этого маленького угла. Пусть а - «новая точка», построенная по сценарию 2. Вторая
точка b смежного ребра (не принадлежащая текущему ребру) сдвигается по формуле:
,
( 10 )
Nа=2 - число ребер, выходящих из точки а; Nb - число ребер, выходящих из точки b;
rа – радиус-вектор точки а; rb - радиус-вектор начального положения точки b; r - радиус-вектор нового положения точки b;
5) Строится новый треугольник, образуемый текущим ребром, одним из граничных рёбер, смежных к нему и имеющий с ним малый угол, и одним новым граничным
ребром, построенным напротив этого маленького угла.
6) Слияние двух смежных граничных рёбер, имеющих общую точку и образующих малый острый угол между собой в одно внутреннее (не граничное) ребро. Две
различающиеся точки смежных ребер (точки а и b) «сливаются» в одну точку. Радиус-вектор получающейся точки считается по следующей формуле:
,
( 11 )
Nа - число ребер, выходящих из точки а; Nb - число ребер, выходящих из точки b;
rа – радиус-вектор точки а; rb - радиус-вектор точки b; r - радиус-вектор образующейся в результате слияния точки;
7) Цикл из четырех ребер разбивается на два треугольника новым построенным пятым ребром, соединяющий противоположные тупые углы четырехугольника.
8) Цикл из четырех ребер разбивается на два треугольника дополнительным пятым ребром, соединяющий «особую» вершину с противоположной ей вершиной.
9) «Особая» вершина цикла из четырех ребер выбрасывается. Из оставшихся трех вершин формируем треугольник.
10)
Построение треугольника закрывающего граничное множество, состоящее из трёх рёбер. Проводиться специальная проверка, чтобы не применять этот метод
для первого шага – когда обрабатывается граничное множество затравочного треугольника.
Рисунок . Определение новой точки поверхности (p4), которая
вместе двумя концевыми точками (p1,p2) ребра (r1) образует
новый треугольник (t2).
Треугольник образуется текущим граничным
ребром и двумя новыми граничными рёбрами.
Пунктиром показаны новые граничные рёбра.
Треугольник образуется текущим и последующим
(предыдущим) граничным рёбрами текущего
граничного набора и одним новым граничным ребром
Пунктиром показано новое граничное ребро
образующее новый треугольник.
Построение треугольника закрывающего
граничное множество, состоящее из трёх рёбер.
Построение треугольника, геометрически
объединяющего или разделяющего
текущие граничные массивы.
Слияние двух граничных рёбер образующих
острый угол между собой в одно внутреннее
ребро.
Определение новой точки поверхности (p4),
которая вместе двумя концевыми точками (p1,p2)
ребра (r1) образует новый треугольник (t2).
Последовательно производя обход граничных
рёбер и добавляя к каждому по новому
треугольнику может возникнуть наложение
треугольников.
Повторяя процедуру обхода граничных рёбер
выполняется послойное построение поверхности.
После каждой процедуры обхода образуется
новый набор граничных рёбер. Граничные рёбра
могут образовывать несколько замкнутых
ломаных.
После обхода всех рёбер первого
треугольника получаем фрагмент
поверхности из 4 треугольников и новое
множество граничных рёбер.
Для случаев наложения треугольников
применяется ряд методов для устранения
конфликтов, например, сшивание двух
треугольников общим ребром.
Построенная поверхность. Поверхность
считается построенной когда не остаётся ни
одного граничного ребра.
Площадь треугольника
1    
Str  ( p2  p1 )  ( p3  p1 )
2
После этого для каждого узла сетки триангуляции производиться
суммирование площадей тех треугольников, для которых
данный узел является вершиной. Полученная суммарная
площадь делиться на три, полученный результат запоминается
как площадь поверхностного элемента с координатами и
нормалью соответствующих данному узлу сетки триангуляции.
Площадь сферического
треугольника
Площадь тороидального
треугольника
O
Y
G
В L
С
D
K
А
E
F
O
X
площадь поверхностного
элемента
После этого для каждого узла сетки
триангуляции производиться суммирование
площадей тех треугольников, для которых
данный узел является вершиной. Полученная
суммарная площадь делиться на три,
полученный результат запоминается как
площадь поверхностного элемента с
координатами и нормалью соответствующих
данному узлу сетки триангуляции.
Входные данные: координаты и радиусы атомов, сферических фрагментов, параметры
тороидальных фрагментов, данные о разбиении пространства на кубические области,
координаты проецируемой точки.
Определение номера
кубической области, в
которой лежит точка.
Определение ближайшего к данной
точке атома из находящихся в
найденном кубе.
Для найденного ближайшего атома – проверка: входит ли
проецируемая точка в один из запрещённых конусов.
Есть
вхождени
е?
Проецирование на сферу
атома – по формулам (1),(2)
Выход с выдачей
координат точки проекции.
Рисунок 40. Блок схема алгоритма проецирования
Проверка: если
ближайшего атома
найдено не было –
поиск среди всех
атомов.
Определение соседних с
текущим атомом
тороидальных фрагментов.
Перебор всех
соседних
тороидальных
фрагментов.
Определение проекции заданной точки на текущий
тороидальный фрагмент.
Проверка угла ` из выражений (3) с интервалом допустимых углов свободной
обкатки () текущего тороидального фрагмента.
Вычисление и запоминание координат,
нормали точки проекции на тороидальный
фрагмент, а также расстояния от исходной
точки до точки проекции.
`()
Обработаны все
тороидальные
фрагменты?
Есть ли хотя бы
одна проекция на
тороидальный
фрагмент.
Выход с выдачей координат точки
проекции ближайшей к исходной
точке.
Цикл по все
вогнутым
сферическим
фрагментам.
Переход к обработке
следующего фрагмента
Определение соседних
вогнутых сферических
фрагментов.
Проверка корректности проецирования на текущую вогнутую сферу – принадлежит
ли проекция сферическому треугольнику, образованному точками контактов сферы
обкатки и опорных атомов.
Вычисление и запоминание
координат, нормали точки проекции
на текущий сферический фрагмент, а
также расстояния от исходной точки
до точки проекции.
Выход с выдачей координат точки проекции
ближайшей к исходной точке.
Есть ли
принадлежность
сферическому
треугольнику?
Обработаны все
вогнутые
сферические
элементы.
Переход к
обработке
следующего
фрагмента.
Начало
SetTris() определение параметров тройных точек первичной
обкатки
CorTris() коррекция тройных точек - удаление совпадающих
SetTwos() определение осей первичной обкатки
CorTrisTw() коррекция данных о тройных точках с учётом
осей.
SetTrNeib() определение масивов тройных точек
окружающих каждый атом
SetTrNbTr()определение тройных точек окружающих каждую
тройную точку
SetTwNeib() определение осей окружения атомов
SetConect() установление массивов поверхностной связности
тройных точек
CorTwos() кррекция данных об осях с учётом поверхностной
сязности(если дополнить процедуру опредления
поверхностной связности то эта функция не нужна)
SetTrsCnf()
заполнение массива конфликтных соседих тройных точек для
каждой тройной точки
SetCnfGrp() формирования набора массов конфликтносязанных тройных точек
SetSecRol() исполнение процедуры вторичной обкатки (с
проверкой и устранением конфликтов)
SetAtomCn() установка набора запрещённых конусов для
каждого поверхностного атома.
SetScForb()
определение параметров, описывающих запрещённые вторичной
обкаткой области
ScStAnlz()
анализ траекторий вторичной обкатки для определения разрывов
связности между парами атомов связаными конфликными торами
первичной обкатки с целью нахождения несвязанных замкнутых
поверхностей
SetCnnAtm()
задания массива поверхнотных атомов для каждой замкнутой
поверхности
Конец
Разбиение полной поверхности на
замкнутые несвязанные
поверхности
• а) Атомы могут быть связанны первичными торами обкатки,
однако при разрыве этих торов связь может сохраняться засчет
торов вторичной обкатки. Строя «дерево» торов первичной и
вторичной обкатки мы находим связанные ими коллективы
атомов. Таким образом, мы находим связанные наборы атомов
(молекулы) и отделяем различные молекулы друг от друга.
• б) Внутренние полости обычно имеют нормаль, направленную
внутрь полости. Написав формулу для нахождения объема
можно видеть, что объем в этом случае получается
отрицательным. Это позволяет отфильтровать такие полости
• с) Внутри полости молекулы может находится другая
«запертая» молекула. Ее поверхность можно отбросит по
следущим признакам – ее внутренние атомы являются также и
внутренними атомами поверхности большой молекулы. При
этом мы выбираем поверхность большего среднего радиуса
Аналитические производные обзор работ
•
•
•
•
Нахождение градиентов параметров поверхностных элементов (эта
задача для алгоритма GEPOL построения поверхности рассмотрена в
[58]
На основе полученных выше градиентов расчитываются градиенты
матриц уравнений и градиенты энергии ( для COSMO [16],[62], для
РСМ [63]) , что для РСМ описанная в [63] методика работает, только
если мы знаем обратную матрицу.
Работа с обратной матрицей создает ряд проблем. Если число
поверхностных элементов велико, то ее размер используемых матриц
велик и не умещается в оперативной памяти. Вызов же из
компьютерной памяти на жестком диске требует огромого времени.
Кроме того обращение матриц [1b], [65] – тяжелая расчетная задача.
Для методов с укрупнеными поверхностными элементами [17]
обратная матрица мала, ее вычисление и использование не
проблематично. Используя методику в [15] (смотр. 3.2.2) можно
итерационно найти обратную матрицу. Используя методы из [63] можно
найти аналитические градиенты
Точки на сфере
Производных от параметров
поверхностных элементов
(координат, нормалей, площадей)
для точек находящихся на
поверхности сферы атома с
координатами центра:


rs  r1
(1)

n s  0
S s  0
- то есть точки перемещаются
вместе с центром атома.
Первичная и вторичная обкатка
двух атомов (сфер).
Торы обкатки.
вторичная обкатка двух сфер
Точки на поверхности
вторичных сфер устойчивых
положений
Первичная и вторичная обкатка для точек
находящихся на сфере шар-зонда при его опоре на
три атома (сферы)


 ( r  r1 ) 2  a 2

 2
 
2
( r  r2 )  b
 ( r  r ) 2  c 2
3


 


(r  r1 )  (r  r1 )  0
 


((
r

r
)


r
)0

2



 ((r  r )  r )  0
3

 

 

(r  r1 )  r  (r  r1 )  r1
 
 
 
 r   [r  r2  r  r3  ]






 
 




r


[
r

r

r

 2
 r3 ]

(r  r1 )  r1


 
 
 

(r  r1 )  [r  r2  r  r3 ]

  

   
(r  r1 )  r1
r        [r  r2  r  r3 ] 
(r  r1 )  [r  r2  r  r3 ]
 

   
(r  r2 )  r2
       [r  r1  r  r3 ] 
(r  r2 )  [r  r1  r  r3 ]
  
   
(r  r3 )  r3
       [r  r2  r  r1 ]
(r  r3 )  [r  r2  r  r1 ]
Площадь треугольника и её
дифференциал.




s  a b

1  1 2
S s 
s
2
2
  11  
1 2 1 1
S   s 
2(s  s ) 
(s  s )
2
4 s 2
4S

   
   
   

   
   
   



a  r2  r1



b  r3  r1
  
   
(s s )  a  b  a  b  a  b  a  b  a  b  a  b  b  a  b  a
a  b a  b 


 a  b  a  b 




 a  a  b  b 

 

 

 a a  b   b a  a  b 
 

 


 a  a b
 a  b b  a 
2
  2
  2
       
1 (a  a )b  (b  b )a  (a  b )(a  b )  (a  b )(a  b )
S 
4
S
 b  a b  a   
 b  b a  b  a a  b 
2
Градиенты площади многоугольных
граничных элементов.
• Градиенты площади граничных
элементов считаются как одна третья от
суммы градиентов площадей
составляющих их треугольников.
Нормали граничных
элементов и их градиенты
N eg

n g ( new ) 

 (a
gi a g ( i 1)

 (a
gi a g ( i 1)
i 1
N eg

sign n g a gi a g( i 1) )

sign n g a gi a g( i 1) )
i 1
Lgi=sign(ng·[agi x ag(i+1)])
sg=∑iLgi[agi x ag(i+1)]
∆jkagi
= ∆jkrgi-∆jkrg
∆jkag(i+1)= ∆jkrg(i+1)-∆ jkrg
Sgi= ([∆jkagi x ag(i+1)]+ [agi x ∆jkag(i+1)]) Lgi
∆jkng(new)=(∑i=1 Neg(Sgi- ng(new) (ng(new)· Sgi)))/|sg|
Неэлектростатическая часть энергии в модели
Абагяна : градиенты
1) поверхностный элемент на первичной сфере, опирающейся на три атома, отображается в точку. Его площадь на SАS
нулевая
2) поверхностный элемент на торе, опирающимся на два атома, отображается в линию. Его площадь на SАS нулевая
3) поверхностный элемент на атоме отображается в сферический тругольник. Его площадь на SАS пропорцианальна его
площади [(ratom+ p_rol)/ ratom]2 Градиент площади по смещению атомов равен нулю
4) граничный поверхностный элемент между тором и первичной сферой, опирающейся на три атома – нулевая площадь и
градиент
5) граничный треугольный поверхностный элемент между тором (и/или первичной сфере, опирающейся на три атома) и
атомом – площадь считается по формуле (187) для трех точек (r1n,r2n,r3n) Градиент площади считаем по формуле (192) , где
градиенты r1n,r2n,r3n:
∆r1n=∆r1+∆n1*p_rol
(1)
∆r2n=∆r2+∆n2*p_rol
(2)
∆r3n=∆r3+∆n3*p_rol
(3)
Для подсчета неэлектростатической части энергии используем формулу Абогяна:
a = r2n-r1n
(6)
b = r3n-r1n
(7)
sj = |axb|/2
(8)
Градиент энергии неэлектростатической части:
∆G=∑j σj ∆sj + b -в ккал/(моль∙А) - j=1,N по всем (с ненулевым градиентом на SAS) поверхностным элементам
S4 = 2|axb|
(11)
∆a
= ∆r2n-∆r1n
(12)
∆b
= ∆r3n-∆r1n
(13)
∆sj
= ((b·b)*(∆a·a)+(a·a)*(∆b·b)-(a·b)(( ∆a·b)+(a·∆b)))/S4 (14)
Проблемы при расчете градиентов и оптимизации с
их использованием
• Вырождение - тройная точка опирается не на три, а
на большее число атомов (ароматические кольца)
• Происходит резкая перестройка поверхности из-за
изменения радиуса вторичной обкатки самой
программой при автоматической его настройке
• Происходит резкая перестройка поверхности из-за
появления или исчезновении тора обкатки
(первичного или вторичного) при сдвиге тройных
точек или изменения узкого перешейка тора на
расстояниях близких к критическому
• Могут быть скачки градиентов или «колебания»
оптимизатора вблизи вышеописанных точек
Проверка, что все атомы лежат внутри одной и
только одной из замкнутых поверхностей.
Оределение внутренних и внешних атомов

 
r  r1  
1
 0 if r1 inside of VS

   3 dS  
4 | r  r1 |

1 if r1 outside of VS
Аналитические производные
для COSMO


1   1
factor 
 in    1

2

Aq   BQ  2 DT Q
q   A1BQ  2 A1DT Q
qT  QT BT A1  2QT DA 1
E  ( factor)QT Dq  ( factor)2QT DA1DT Q  
factor T
q Aq
2
1 dE
d (QT DA1DT Q)
dD
1
dA
dDT
 2
 QT
(2 A1DT Q)  ( 2QT DA1 ) (2 A1DT Q)  (2QT DA1 )
Q
factor dR
dR
dR
2
dR
dR
 QT
dD
1
dA
dD T
dD T
1
dA
q  qT
q  qT
Q  2q T
Q  qT
q
dR
2
dR
dR
dR
2
dR






Матрица D и ее производная

1
1
dij    
2 R r
j
 i





   
1 ( Ri  r j )( Ri  r j )
d ij  
  3
2
Ri  r j
Матрица A и ее производная
i j

1

aij    
 ri  r j






   
(ri  rj )(ri rj)
aij  
  3
ri  rj
i j
aii  Const
1
Si
aii  Const
Si
2 Si3
Аналитические производные
для РСМ
Aq DQ
Aq  BQ

T
q  A BQ
T
q  (A ) D Q
1
1
*
T
T
E  Q Dq  Q DA BQ  q  BQ  Q B q  q  Aq
T
1
T
T

T
T


1
T
dE d (Q DA BQ)
dD
dA
dB

Q
A BQ  (Q DA ) ( A BQ)  (Q DA ) Q 
dR
dR
dR
dR
dR
dD
dA
dB
Q
q  (q )
q  (q )
Q
dR
dR
dR
T
T
T

T
1

T
T
1
1
T
1
Условие нормализации
Матрица D и ее производная
1  1
 
d 

2 R  r




ij
i
j
(1)
   
1 ( R  r )( R  r )
d  
 
2
R r
i
j
i
j
3
ij
i
j
Матрица B и ее производная
  
n r  R 
b    S * ( Feps )
r R
i
i
j
3
ij
i
i
j

  
  
  
 n r  R 







r  R
n
r

R
n
r

R
n
b  ( Feps)
 S 
 S 
 S 3

 r R
r

R
r

R
r

R

i
i
j
3
ij
i
j
i
i
j
3
i
i
i
i
j
3
i
j
4π ∙FEPS=2(1- ε)/(ε+1)
i
j
i
i
i
i
j
j
5
r  R r R  S 
i
j
i
j
i


Матрица A и ее производная
i j
  
n r  r 
a    S * ( Feps)
r r
i
i
j
3
ij
i
i
j
  
  
  
 n r  r 
nr  r 
n r  r 
n r  r      
a  ( Feps)   S     S    S  3   r  r r  r S 
r r
r r
r r
 r r

i
i
j
3
ij
i
i j
i
i
j
i
3
i
j
i
i
j
3
i
j
i
j
i
i
j
5
i
i
i
j
i
j
i
j
2
 1
Dcc5=
  
n r  r 
a  Dcc5     S * ( Feps)  Dcc5   a
r r
i
ii
i
j
3
i j
i
a   a
ii
i j
ij
j
i
i j
ij
Дифференцирование SGB
1 1
G 
2  in

1  1
 


 
 i, j

Qi  Q j
1

2
 2
Ri , j
Ri , j  ai  a j  exp( c a a
i
gi , j 
1
j
1 1
 
1
2  in
1
2
)
G  

2
 2
Ri , j
Ri, j  ai  a j  exp( c  )
i
1 1

2  in

1  1
 


1  1
 


1
 Q Q  g
 i, j i j i, j
1

1
2

1
   Q  Q  g
i
j
i, j
 i, j
1
1

2
j






1

g i , j   
 2
  2

R
 R i, j  a i  a j  exp( i , j ) 
c a i a j



 Ri, j

 2

  1
R i, j
 g i , j3  



2( R i, j  R i, j ) 1  exp( c a a )   a i a j  a i a j
 c
i j 
2 














   Ri, j  Ri, j  Ri, j  Ri, j  Ri, j  Ri, j  2( Ri, j  Ri, j )

2
 2
 2 

R i, j 
R i, j
1 
 exp(
)
c

a

a
c a i a j 

i j 




Градиенты борновских радиусов
ai 
   

 (ns  (rs  Ri ))dS 
i
I n  

  n


rs  Ri


1


2  An  I ni  A0 
n

ai 
  An  I ni
2
 n

2  An  I ni  A0 
n

n 3
   

 (ns  (rs  Ri )) S s 
i
I n  

  n
s

rs  Ri


n4
1
n3
I ni

 
J ni i

Jn
n 3
4 n
n3

J ni
1

n  3 I i n 4
n
J si ,n

n3
n4
 


J ni   J si ,n



 
 


  

  
  
  

 ns  (rs  Ri )  ns  (rs  Ri ) S s  ns  (rs  Ri ) S s n ns  (rs  Ri ) (rs  Ri )  (rs  Ri ) S s 



  n
  n2


rs  Ri
rs  Ri



1
n4
s
  

    n        n / 2  n       ( n  2) / 2 
 rs  Ri    (rs  Ri )(rs  Ri )
     (rs  Ri )(rs  Ri )
  2 (rs  Ri )(rs  Ri ) 
 2 


 
  

n (rs  Ri )  (rs  Ri )

  n 2
rs  Ri

 
 J ni
  
(
n
 (r  R )) S
J si ,n  s s  in s

rs  Ri
  
(
n
 (r  R )) S
J ni   s s  in s   J si , n

s
s
rs  Ri
1
1




Energy PCM vs COSMO as a function of the grid step (for ε=∞)
Y axis – |[E(cosmo)-E(pcm)]/E(pcm)|
0,12
0,1
0,08
0,06
Series
1
0,04
0,02
0
0
0,1
0,2
0,3
0,4
0,5
Gradiens PCM vs COSMO (for ε=∞)
Gradiens COSMO and SGB vs PCM (for ε=78.5)
относительная ошибка
Относительная ошибка аналитических градиентов COSMO и SGB в
сравнении с PCM
0,4
0,375
0,35
0,325
0,3
0,275
0,25
0,225
0,2
0,175
0,15
0,125
0,1
0,075
0,05
0,025
0
COSMO
SGB
0
2
4
6
8
10
12
абсолютная величина аналитического градиента, PCM
14
16
Gradiens PCM, COSMO, SGB: anal vs num (for ε=78.5)
относительная ошибка
Относительная ошибка численного и аналитического
градиента
0,4
0,375
0,35
0,325
0,3
0,275
0,25
0,225
0,2
0,175
0,15
0,125
0,1
0,075
0,05
0,025
0
COSMO
PCM
SGB
1
2
3
4
5
6
7
8
9 10 11 12 13 14 15 16
абсолютная величина аналитического градиента,
(ккал/моль А)