分子动力学后处理自编程系列（4）---氢键统计

氢键是一种特殊的静电作用，是除范德华力外的另一种分子间作用力，在小分子(水、离子溶液等)及聚合物(DNA、蛋白质、无机材料等) 的三维结构和特性方面起着重要作用。氢键的大小介于化学键与范德华力间，不属于化学键，但有键长、键能，氢键具有饱和性、方向性。本程序利用LAMMPS输出的dump文件，统计水中氢键数量。

本程序的优势主要表现在：
💖输入明确：dump文件（包含Type、mol-ID以及x，y，z）
💖简单易修改：程序思路清晰，可根据自己需求在本程序基础上进行修改实现特性需求计算，例如氢键键长分布、键角分布等。

1、编程思路

（1）读取dump文件，将每一帧的原子信息存储为变量atom_data；
（2）在帧循环条件下，按照O和H两种原子将data信息拆分；
（3）选出潜在的能够当做供体的O原子，并确定与氧原子相连接的H原子；
（4）根据O-O之间的距离大于2.5，小于R判断供体O对应的潜在的能够当做受体的O原子；
（5）循环供体O，计算该供体O连接的H与潜在受体O之间的夹角α 来判断是否形成氢键；
（6）记录所有帧下氢键数目并写出供体、受体信息；
（7）根据需要完成特性信息筛选。

2、计算步骤

(1) 读取dump文件

dump的读取和该系列之前的读取代码一致，请移步之前的教程，复制该部分代码，为了避免教程冗长，在后续的教程中该部分一并省略。

(2) 几何判断参数设置

% 输入判断条件
R = 3.5;                   	% 几何规则距离判断条件
aplha = 30;                 % 几何规则角度判断条件

(3) 选出形成氢键相关的O和H原子

% Type 1 ------ O
% Type 2 ------ H

for i = 1 : frame_num
    n = 0;                                          
    m = 0;                                          
    for j = 1 : num_atoms
        if (atom_data(j,2,i) == 1)                  
            n = n + 1;
            data_Oxygen(n,:,i) = atom_data(j,:,i);
            XYZ_Oxygen(n,:,i) = atom_data(j,3:5,i);
            Type_Oxygen(n,:,i) = atom_data(j,2,i);
        elseif (atom_data(j,2,i) == 2)
            m = m + 1;
            data_Hydrogen(m,:,i) = atom_data(j,:,i);
            XYZ_Hydrogen(m,:,i)  = atom_data(j,3:5,i);                         
            Type_Hydrogen(m,:,i) = atom_data(j,2,i);                           
        end
    end
end


% No.3 选出每个donor对应的H
for i = 1 : frame_num                                                       
    for j = 1 : n                                                                                                         
        k = 0;                                                              
        for v = 1 : m                                                          
        	dx = XYZ_Oxygen(j,1,i) - XYZ_Hydrogen(v,1,i);
        	dy = XYZ_Oxygen(j,2,i) - XYZ_Hydrogen(v,2,i);
            dz = XYZ_Oxygen(j,3,i) - XYZ_Hydrogen(v,3,i);
            distance = sqrt(dx^2 + dy^2 + dz^2);
            if ( 0.90<distance && distance< 1.10)
                k = k + 1;
                donor_H(k,:,j,i) = data_Hydrogen(v,:,i);                         
            end
        end
    end
end

(4) 选出每个donor的acceptor

% No.4 选出每个donor的acceptor
for i = 1 : frame_num
    for j = 1 : n                                                           
       y = 0;                                                               
       for x = 1 : n                                                        
          d_x = XYZ_Oxygen(j,1,i) - XYZ_Oxygen(x,1,i);
          d_y = XYZ_Oxygen(j,2,i) - XYZ_Oxygen(x,2,i);
          d_z = XYZ_Oxygen(j,3,i) - XYZ_Oxygen(x,3,i);
          distance_0 = sqrt(d_x^2 + d_y^2 + d_z^2);
          if ( distance_0 < R )  
              y = y + 1;
              donor_A(y,:,j,i) = data_Oxygen(x,:,i);
          end         
       end
   end
end

(5) 氢键数目统计

% No.5 氢键初统计
for i = 1 : frame_num
    n_hbonds(i) = 0;											% 氢键数
    for j = 1 : n                                                           
        XYZ_donor_now = XYZ_Oxygen(j,:,i);
        n_H(j,i) = sum(donor_H(:,1,j,i)~=0);                               
        for ii = 1 : n_H(j,i)                                               
            XYZ_H_now = donor_H(ii,3:5,j,i);
            O_H = XYZ_H_now - XYZ_donor_now;
            n_A(j,i) = sum(donor_A(:,1,j,i)~=0);                            
            for jj = 1 : n_A(j,i)                                           
                XYZ_A_now = donor_A(jj,3:5,j,i);
                O_O = XYZ_A_now - XYZ_donor_now;
                sita = (180/pi)*acos(dot(O_H,O_O)/(norm(O_H)*norm(O_O)));

                if (sita < aplha)
                    n_hbonds(i) = n_hbonds(i) + 1;
                    % 存储氢键的信息，行数为氢键数，
                    % 第一列为供体TYPE，第二列为供体H的TYPE，第三列为受体TYPE，第四列为供体原子ID，第五列为供体H原子ID，第六列为受体原子ID；
                    Hbonds_info(n_hbonds(i),1,i) = Type_Oxygen(j,1,i);
                    Hbonds_info(n_hbonds(i),2,i) = donor_H(ii,2,j,i);
                    Hbonds_info(n_hbonds(i),3,i) = donor_A(jj,2,j,i);                    
                                               
                    Hbonds_info(n_hbonds(i),4,i) = data_Oxygen(j,1,i);
                    Hbonds_info(n_hbonds(i),5,i) = donor_H(ii,1,j,i);
                    Hbonds_info(n_hbonds(i),6,i) = donor_A(jj,1,j,i);               
                    
                end
            end
        end
    end
end

(6) 结果验证

3、结语

后续我们会开发一个免费的氢键统计APP，功能强大，操作简单，开发完成会在本账号和微信公众号“分子模拟CLUB”上发布，敬请期待。