DFTBで分子のオンサイトエネルギーを表示する簡単なコード

本コードはハミルトニアンを出力してることが前提となります。1番目の引数に分子座標のファイル。2番目の引数を0にしたらs軌道、1ならpx軌道、2ならpy軌道みたいな感じです。

改良したら、指定した原子のオンサイトエネルギーを出したりと色々できそう。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#define NN 4000

int main (int argc, const char * argv[]) {
__    if(argc == 1){
__    printf(“Usage: \n[your_name home] dftb-site mol.xyz <0,1,2,3>  \n”);
__    printf(” \n”);exit(0);
}
FILE *fin; // mount point
int atom_num; // the number of atoms
char element[NN][10];
double x[NN];
double y[NN];
double z[NN];
double on_site[NN];
double site;
char str[256];
char path[256];
int  element_num[NN];
int  chosen_element_num[NN];
char chosen_element[NN][10];
int  atom_index[NN];
int  basis_per_atom[NN];
int  start_index[NN];
int  show_list[NN];

int  c[120];
int  flag_static=0;
int  flag_s=0;
int  flag_px=0;
int  flag_py=0;
int  flag_pz=0;

int  NBasis,i,j,k,i1,i2,op,cnt,sum;

if( (fin=fopen(argv[1],”r”)) == NULL ){  // File opening
printf(“Can’t open file!\n”);
exit(1);
}
strcpy(str,argv[2]);
if(strcmp(str,”0″)==0){

flag_s=1; op=0;

} else if(strcmp(str,”1″)==0){

flag_px=1;op=1;

} else if(strcmp(str,”2″)==0){

flag_py=1; op=2;

} else if(strcmp(str,”3″) ==0){

flag_pz=1; op=3;

}

if( fscanf(fin,”%d”,&atom_num) != EOF){
}else{
printf(“This file may not be xyz format\n”); exit(1);
}

for(i=0;i<atom_num; i++){
if( fscanf(fin,”%s”,element[i]) == EOF ){ printf(” Something wrong for element reading \n”); }
if( fscanf(fin,”%lf”,&x[i]) == EOF ){ printf(“Something wrong for %d-th x reading\n”,i); }
if( fscanf(fin,”%lf”,&y[i]) == EOF ){ printf(“Something wrong for %d-th y reading\n”,i); }
if( fscanf(fin,”%lf”,&z[i]) == EOF ){ printf(“Something wrong for %d-th z reading\n”,i); }
}
fclose(fin);

for(i=0;i<atom_num;i++){    show_list[i]=0;    }
cnt=0; sum=0;
for(i=0;i<atom_num; i++){
if(strcmp(element[i],”C”) == 0 )      {  basis_per_atom[i]=4; show_list[cnt] = sum +op; cnt++; sum += 4; }
else if(strcmp(element[i],”N”) == 0 ) {  basis_per_atom[i]=4; show_list[cnt] = sum +op; cnt++; sum += 4; }
else if(strcmp(element[i],”O”) == 0 ) {  basis_per_atom[i]=4; show_list[cnt] = sum +op; cnt++; sum += 4; }
else if(strcmp(element[i],”H”) == 0 ) {  basis_per_atom[i]=1; show_list[cnt] = sum   ;  cnt++; sum += 1; }
}

if( (fin=fopen(“hamsqr1.dat”,”r”)) == NULL ){  // File opening
printf(“Can’t open hamsqr1.dat!\n”);
exit(1);
}

fgets(str,256,fin);
if( fscanf(fin,”%s”,str) == EOF ){ printf(“error at 83\n”,i); }
if( fscanf(fin,”%d”,&NBasis) == EOF ){ printf(“error at 84\n”,i); }
fgets(str,256,fin);

if( fscanf(fin,”%s”,str) == EOF ){ printf(“error at 83\n”,i); }

fgets(str,256,fin);
fgets(str,256,fin);
fgets(str,256,fin);
printf(“NBasis=%d\n”,NBasis);

for(i=0;i<NBasis;i++){
for(j=0;j<NBasis;j++){
if( fscanf(fin,”%lf”,&site) == EOF ){ printf(“Something wrong for %d-th z reading\n”,i); }
if( i == j ){
on_site[i]=27.21134*site;
printf(“on_site=%lf\n”,on_site[i]);
}
}
}

// results

printf(“%d\n\n”,atom_num);
for(i=0;i<atom_num;i++){
j=show_list[i];
printf(“%s  %lf  %lf  %lf  %lf\n”,element[i],x[i],y[i],z[i],on_site[j]);
}
return 0;
}

poissondd について

日々感じたことを書いていこうと思います。 テーマは分子エレクトロニクスから、MemristorのBehavioral modelまでいろいろ。
カテゴリー: 未分類 パーマリンク

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト /  変更 )

Google フォト

Google アカウントを使ってコメントしています。 ログアウト /  変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト /  変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト /  変更 )

%s と連携中