| 查看: 1006 | 回复: 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 |
回复此楼» 猜你喜欢
面上再次挂了，太难了，躺也躺不了，倦也卷不过，小学校之殇！
已经有25人回复
-
2026年国自然面上资助率
已经有22人回复
-
面上提前没消息，有中的吗
已经有15人回复
-
基础研究怎么拉横向，学校到款任务越来越多，难以完成 拉横向，都有哪些途径啊
已经有12人回复
-
这种情况还有戏吗
已经有12人回复
-
微信指数没变化，科研之友没阅读
已经有19人回复
-
准备明年的基金了
已经有7人回复
-
2027年申博
已经有3人回复
-
一个有机合成实验室都需要哪些设备？
已经有7人回复
-
HXDI做水性聚氨酯乳液，是不是特别容易出渣
已经有3人回复
2楼2013-05-08 17:09:59