| 查看: 840 | 回复: 1 | |||
Joannaouc木虫 (著名写手)
|
[交流]
利用Grimme D2 scheme 手动计算dispersion (by MATLAB)已有1人参与
|
|
好久之前写的东西，可能对大家有用。欢迎自行修改。请轻拍砖。 P.S. Dmol5.0 (not quite sure) 以上和VASP5.2以上可以利用程序在结构驰豫的时候就把dispersion加进去。这个当初是针对没有dispersion的dmol写的，刚写完没多久程序就更新到自带了。 function [Edisp] = Edisperse(filename) %define the function to calculate dispersion energy %modified on 09-03-2011 %open file and read fid = fopen(filename,'r'); %following is the matrix storing Ro and C6 information %of the atoms in the system % %first column is C6, second column is Ro % %line 1 for H para(1,1) = 0.14; para(1,2) = 1.001; %line 2 for C para(2,1) = 1.75; para(2,2) = 1.452; %line 3 for N para(3,1) = 1.23; para(3,2) = 1.397; %line 4 for Li para(4,1) = 1.61; para(4,2) = 0.825; %line 5 for Ti para(5,1) = 10.80; para(5,2) = 1.562; %following is some other parameters s6 = 0.75; d = 20; %read in coordinates and atom type of the atoms [atype, x, y, z] = textread(filename, '%s %f %f %f'); l=length(x); %allocate integer number to identify each atom kind for i=1:l if (strcmp(atype(i),'H')==1) atomtype(i) = 1; elseif (strcmp(atype(i),'C')==1) atomtype(i) = 2; elseif (strcmp(atype(i),'N')==1) atomtype(i) = 3; elseif (strcmp(atype(i),'Li')==1) atomtype(i) = 4; else atomtype(i) = 5; end end atomtype = atomtype'; Edisp = 0; for i = 1: (l-1) for j = (i+1):l Rr = para(atomtype(i),2) + para(atomtype(j),2); C6ij = sqrt (para(atomtype(i),1).* para(atomtype(j),1)); Rij = sqrt ((x(i)-x(j)).^2+(y(i)-y(j)).^2+(z(i)-z(j)).^2); fdmpij = 1./(1+exp(-d*((Rij./Rr)-1))); Edisp = Edisp + (C6ij .* fdmpij) ./ (Rij.^6); end end Edisp = - Edisp .* s6; %unit here is k kJ/mol Edisp = Edisp .* 1000; %unit here is kJ/mol Edisp = Edisp ./ 2625.5 .*27.211; %unit here is eV |
回复此楼» 猜你喜欢
投稿Elsevier的Neoplasia杂志，到最后选publishing options时页面空白，不能完成投稿
已经有19人回复
-
职称评审没过，求安慰
已经有19人回复
-
垃圾破二本职称评审标准
已经有12人回复
-
EST投稿状态问题
已经有7人回复
-
谈谈两天一夜的“延安行”
已经有15人回复
-
毕业后当辅导员了，天天各种学生超烦
已经有4人回复
-
聘U V热熔胶研究人员
已经有10人回复
-
求助文献
已经有3人回复
-
投稿返修后收到这样的回复，还有希望吗
已经有8人回复
-
三无产品还有机会吗
已经有6人回复
2楼2013-05-08 17:09:59