| 查看: 845 | 回复: 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 |
回复此楼» 猜你喜欢
全日制（定向）博士
已经有5人回复
-
假如你的研究生提出不合理要求
已经有10人回复
-
萌生出自己或许不适合搞科研的想法，现在跑or等等看？
已经有4人回复
-
Materials Today Chemistry审稿周期
已经有4人回复
-
参与限项
已经有3人回复
-
实验室接单子
已经有4人回复
-
对氯苯硼酸纯化
已经有3人回复
-
求助:我三月中下旬出站，青基依托单位怎么办？
已经有12人回复
-
所感
已经有4人回复
-
要不要辞职读博？
已经有7人回复
2楼2013-05-08 17:09:59