| 查看: 1017 | 回复: 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 |
回复此楼» 猜你喜欢
奇怪，两个人的filecode固定段从头到尾一模一样
已经有5人回复
-
帮忙看看fileCode
已经有11人回复
-
FileCode能看出啥？
已经有16人回复
-
小木虫上这么多卖论文的，真有人买论文么？感觉没必要啊
已经有4人回复
-
为什么网上很多人说本周 12号出结果
已经有4人回复
-
基金中了
已经有15人回复
-
filecode
已经有8人回复
-
确定了，国自然21号放榜
已经有7人回复
-
综述论文作为代表作会不会影响评审专家的印象分？
已经有11人回复
-
关于代码变化问题，想知道的进来
已经有24人回复
2楼2013-05-08 17:09:59
100