豆粑粑 matlab不同的方式画分子的表面。

meatball1982

是让孩子们作的，我给整一下先。
这两天写本子，画几个小图，放松一下。

你的有一个分子。
我是从pdb 里把位置给提出来。
设置一个球，以原子位置为球心，放球。
在空间中有很多的点。

通过boundary(真是个好函数，比较新的版本才有）。实现2D和3D边界的显示。

1. clear all
   
2. clc
   
3. clf
   
4. 
   
5. 
   
6. 
   
7. dat=load('../file_data/pos_ala.txt');
   
8. 
   
9. x=dat(:,1);
   
10. y=dat(:,2);
    
11. z=dat(:,3);
    
12. 
    
13. [sx,sy,sz]=sphere(50);
    
14. 
    
15. r=2.2;
    
16. 
    
17. 
    
18. 
    
19. n_p=length(x);
    
20. 
    
21. x_all=[];
    
22. y_all=[];
    
23. z_all=[];
    
24. 
    
25. 
    
26. for i =1:n_p
    
27. %     rat = 1+0.03*i;
    
28.     rat = 1;
    
29.     sx_new=rat*r.*sx;
    
30.     sy_new=rat*r.*sy;
    
31.     sz_new=rat*r.*sz;
    
32.     x_all=[x_all;dat(i,1)+sx_new(:)];
    
33.     y_all=[y_all;dat(i,2)+sy_new(:)];
    
34.     z_all=[z_all;dat(i,3)+sz_new(:)];
    
35. end
    
36. 
    
37. %
    
38. 
    
39. axis equal
    
40. 
    
41. subplot(2,2,1)
    
42. plot3(x_all,y_all,z_all,'.','markersize',3);
    
43. view(90,17.6)
    
44. axis off
    
45. axis tight
    
46. 
    
47. subplot(2,2,2)
    
48. [k_all3 ]=boundary([x_all,y_all,z_all],0.8);
    
49. 
    
50. trisurf(k_all3,x_all,y_all,z_all,'facecolor','red','facealpha',0.2,'edgecolor','y','edgealpha',0.2)
    
51. view(90,17.6)
    
52. axis off
    
53. axis tight
    
54. 
    
55. subplot(2,2,3)
    
56. lin = linspace(-4,4,8);
    
57. lin =[  -3 -2.5 -2 -1.5 -1 -0.5 0 1 2 3 4 5 6 7 8 9 10 ]
    
58. n_l =  length(lin)
    
59. 
    
60. 
    
61. col_mm = jet(n_l);
    
62. hold on
    
63. for i = 1:n_l
    
64.     lb = lin(i)-0.2;
    
65.     ub = lin(i)+0.2;
    
66.     ind = z_all > lb & z_all < ub;
    
67.     tm_x = x_all(ind);
    
68.     tm_y = y_all(ind);
    
69.    
    
70.     k_all=boundary([tm_x,tm_y],0.99-0.02*i);
    
71.    
    
72.    
    
73.     plot3(tm_x(k_all),tm_y(k_all),lin(i)*ones(size(k_all)),'-','linewidth',1+0.3*i,'color',col_mm(i,:))
    
74.    
    
75. end
    
76. 
    
77. 
    
78. view(90,17.6)
    
79. axis equal
    
80. axis off
    
81. axis tight
    
82. 
    
83. subplot(2,2,4)
    
84. lin = linspace(-4,4,8);
    
85. lin =[  -3 -2.5 -2 -1.5 -1 -0.5 0 1 2 3 4 5 6 7 8 9 10 ]
    
86. n_l =  length(lin)
    
87. 
    
88. 
    
89. col_mm = jet(n_l);
    
90. hold on
    
91. for i = 1:n_l
    
92.     lb = lin(i)-0.2;
    
93.     ub = lin(i)+0.2;
    
94.     ind = z_all > lb & z_all < ub;
    
95.     tm_x = x_all(ind);
    
96.     tm_y = y_all(ind);
    
97.    
    
98.     k_all=boundary([tm_x,tm_y],0.99-0.02*i);
    
99.    
    
100.    
     
101.     h=patch(tm_x(k_all),tm_y(k_all),lin(i)*ones(size(k_all)),col_mm(i,:),'edgecolor','none')
     
102.     alpha(h,0.8-0.02*i);
     
103.    
     
104.    
     
105. end
     
106. 
     
107. 
     
108. view(90,17.6)
     
109. axis equal
     
110. axis off
     
111. axis tight
     
112. 
     
113. 
     
114. h=gcf;
     
115. fi_na = '../file_imgs/fig_cont_ala_4sub';
     
116. fun_work_li_035_myfig_out(h,fi_na,3)
     

复制代码

meatball1982

思路是通过每个原子的vdw半径向外进行扩展，在球面，生成均匀颁的散点。
对于不同半径的球面，分别生成散点。

其中，第i个原子的球面的散点要满足距离atom_i 相比于距离atom_j要近。否则，放弃散点。

球面上生成均匀散点是根据文献实现。
E. B. Saff and A. B. J. Kuijlaars, "Distributing Many Points on a Sphere", The Mathematical Intelligencer, 19, 5-11 (1997)

生成网格点是根据tinker中potential.f中的子函数
potgrid
实现。