岐路

日本ではアカデミアでやれるだけの要領の良さはないと言われ、
かといってドイツではテニュアはabsolute zeroと言われたりもしたのだが、
結果的には今の選択は間違っていなかったと思う。

こうなるよりほかなかったのだ。終わったことをどうこう言ったって
全く生産的ではないし、未来が変わるわけでもない。
選んだ道を信じ(信じて選んだ道だが…)、前を見据え進むだけだ。

カテゴリー: 未分類 | コメントをどうぞ

ルーマニアの夢

数日前、ルーマニアの陰鬱な街を歩いてて、
前から来た若い男二人組にすれ違いざまに
手に持ってたチーズ入りパンの最後の一口を
強引に奪われて食べられたけど、
強盗にも会わずそれっきりでした。

という夢を見た

何を示唆しているのだろう。。。?
ルーマニアは、昨日モルドバの歴史について調べていたからそれつながりだろう。
チーズは昨日食べたからそのせいだ
ルーマニア=強盗というイメージで、チーズ入りパンを奪われたんだろう。
しかしみみっちい強盗だ

カテゴリー: 未分類 | コメントをどうぞ

Gromacs(ver. 5.1.4)をMacBookpro にインストールしてみた

条件

OS X El Capitan (Version 10.11.6)
MacBook Pro
Processor 2.6 GHz Intel Core i7
Memory 8GB 1600 MHz DDR3

http://manual.gromacs.org/documentation/5.1.4/install-guide/index.html に従いインストール。

(1) 本家よりgromacs-5.1.4.tarをダウンロード

(2) ~/workspace/ で展開した。よってgrimacesのホームは~/workspace/gromacs-5.1.4/

(3) cd workspace/gromacs-5.1.4

(4) mkdir build

(5) cd build

(6) cmake ..  -DGMX_BUILD_OWN_FFTW=ON -DREGRESSIONTEST_DOWNLOAD=ON -DGMX_MPI=on

— Build files have been written to: /Users/my_home/workspace/gromacs-5.1.4/build と出た。

(7) make

(8) make check

(9) sudo make install

(10) source /usr/local/gromacs/bin/GMXRC

できた。

追記:

which gmx_mpi でgmx_mpiへのソースというかパスがどこにあるか確認。
/usr/local/gromacs/bin/gmx_mpi
Macportsでもgmx_mpiをなぜか入れてて、新しくターミナルを開けるとそっちにパスが言ってしまうので、次の1行をホームの.bashrcに加えた。

echo  “source /usr/local/gromacs/bin/GMXRC” >> ~/.bashrc
そして更新した。
source ~/.bashrc

追記2:
ラボの据え置きMacではなぜかmake checkでこけた。よくわからんでいじくりまわして、
どうやらMPICHの環境がcompatibleでないっぽい。そこで、Macportsでopenmpiを入れて、
sudo port install openmpi-gcc5

mpiがopenmpi-gcc5-fortran意味することをシステムに認知させ(←業界では正しくはどういう表現をするのかがわからなくてぎこちなくなった…)
sudo port select –set mpi openmpi-gcc5-fortran

さらに~/.bashrcに以下のような行を書き加えて、
export DYLD_LIBRARY_PATH=/opt/local/lib/openmpi-gcc5
export LD_LIBRARY_PATH=/opt/local/lib/openmpi-gcc5

export
LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/local/lib:/usr/lib/:/usr/local/lib:/usr/local/bin

上記の(4)-(8)をやり直したらうまく言った。

カテゴリー: 未分類 | コメントをどうぞ

Gaussianのstandard orientationをlogファイルの底の行から抜き出すスクリプト。

#! /usr/bin/bash

echo $#

# this code searches standard orientation from the bottom line of file and prints out the cartecian coordinate.

if [ $# -eq 0 ]; then

     echo   Usage: print-standard mol.log
     echo this code searches standard orientation from the bottom line  of file and print the cartecian coordinate.

     exit

fi

str=$(tac $1 |  grep -n "Standard orientation:" $1 |  tail -n 1 |  awk '{print $1}' )

#str=$(grep -n "Standard orientation:" $1   tail -n 1   awk '{print $1}' )

start=$(echo "${str//:}")

end=$(wc -l $1 |  awk '{print $1}')

#echo $end

end=$(sed ''"$start"','"$end"'!d' $1 |  nl -ba -v$start |  grep "\-\-\-\-\-" |   head -n 3 | tail -n 1  | awk '{print $1}')

#echo $end

N=$(echo $end-$start-5|  bc)

echo $N

echo " " 

#sed ''"$(echo $start+5  bc)"','"$(echo $end-1 bc)"'!d' $1   awk '{print $2,$4,$5,$6}' | sed -e 's/^1 / H /g' -e 's/^6 / C /g' -e 's/^7 / N /g' -e 's/^8 / O /g' -e 's/^16 / S /g' -e 's/^26 / Fe /g' 

#sed ''"$(echo $start+5 | bc)"','"$(echo $end-1| bc)"'!d' $1 |  awk '{print $2,$4,$5,$6}' 

sed ''"$(echo $start+5 | bc)"','"$(echo $end-1| bc)"'!d' $1 |  awk '{print $2,$4,$5,$6}' | sed -e   's/^1 /  H  /g'  -e   's/^2 /  He  /g'   -e   's/^3 /  Li  /g'   -e   's/^4 /  m  /g'   -e   's/^5 /  B  /g'   -e   's/^6 /  C  /g'   -e   's/^7 /  N  /g'   -e   's/^8 /  O  /g'   -e   's/^9 /  F  /g'   -e   's/^10 /  Ne  /g'   -e   's/^11 /  Na  /g'   -e   's/^12 /  Mg  /g'   -e   's/^13 /  Al  /g'   -e   's/^14 /  Si  /g'   -e   's/^15 /  us  /g'   -e   's/^16 /  S  /g'   -e   's/^17 /  Cl  /g'   -e   's/^18 /  Ar  /g'   -e   's/^19 /  K  /g'   -e   's/^20 /  Ca  /g'   -e   's/^21 /  Sc  /g'   -e   's/^22 /  Ti  /g'   -e   's/^23 /  V  /g'   -e   's/^24 /  Cr  /g'   -e   's/^25 /  Mn  /g'   -e   's/^26 /  Fe  /g'   -e   's/^27 /  Co  /g'   -e   's/^28 /  Ni  /g'   -e   's/^29 /  Cu  /g'   -e   's/^30 /  Zn  /g'   -e   's/^31 /  Ga  /g'   -e   's/^32 /  Ge  /g'   -e   's/^33 /  As  /g'   -e   's/^34 /  Se  /g'   -e   's/^35 /  Br  /g'   -e   's/^36 /  Kr  /g'   -e   's/^37 /  Rb  /g'   -e   's/^38 /  Sr  /g'   -e   's/^39 /  Y  /g'   -e   's/^40 /  Zr  /g'   -e   's/^41 /  Nb  /g'   -e   's/^42 /  Mo  /g'   -e   's/^43 /  Tc  /g'   -e   's/^44 /  Ru  /g'   -e   's/^45 /  Rh  /g'   -e   's/^46 /  Pd  /g'   -e   's/^47 /  Ag  /g'   -e   's/^48 /  Cd  /g'  -e   's/^49 /  In  /g'  -e   's/^50 /  Sn  /g'  -e   's/^51 /  Sb  /g'  -e   's/^52 /  Te  /g'  -e   's/^53 /  I  /g'  -e   's/^54 /  Xe  /g'  -e   's/^55 /  Cs  /g'  -e   's/^56 /  Ba  /g'  -e   's/^57 /  La  /g'  -e   's/^58 /  Ce  /g'  -e   's/^59 /  Pr  /g'  -e   's/^60 /  Nd  /g'  -e   's/^61 /  Pm  /g'  -e   's/^62 /  Sm  /g'  -e   's/^63 /  Eu  /g'  -e   's/^64 /  um  /g'  -e   's/^65 /  Tb  /g'  -e   's/^66 /  Dy  /g'  -e   's/^67 /  Ho  /g'  -e   's/^68 /  Er  /g'  -e   's/^69 /  Tm  /g'  -e   's/^70 /  Yb  /g'  -e   's/^71 /  Lu  /g'  -e   's/^72 /  Hf  /g'  -e   's/^73 /  Ta  /g'  -e   's/^74 /  W  /g'  -e   's/^75 /  Re  /g'  -e   's/^76 /  Os  /g'  -e   's/^77 /  Ir  /g'  -e   's/^78 /  Pt  /g'  -e   's/^79 /  Au  /g'  -e   's/^80 /  Hg  /g'  -e   's/^81 /  Tl  /g'  -e   's/^82 /  Pb  /g'  -e   's/^83 /  Bi  /g'  -e   's/^84 /  Po  /g'  -e   's/^85 /  At  /g'  -e   's/^86 /  Rn  /g'  -e   's/^87 /  Fr  /g'  -e   's/^88 /  Ra  /g'  -e   's/^89 /  Ac  /g'  -e   's/^90 /  Th  /g'  -e   's/^91 /  Pa  /g'  -e   's/^92 /  U  /g'  -e   's/^93 /  Np  /g'  -e   's/^94 /  Pu  /g'  -e   's/^95 /  Am  /g'  -e   's/^96 /  Cm  /g' -e   's/^97 /  Bk  /g'  -e   's/^98 /  Cf  /g'  -e   's/^99 /  Es  /g'  -e   's/^100 /  Fm  /g'  -e   's/^101 /  Md  /g'  -e   's/^102 /  No  /g'  -e   's/^103 /  Lr  /g'  -e   's/^104 /  Rf  /g'  -e   's/^105 /  Db  /g'  -e   's/^106 /  Sg  /g'  -e   's/^107 /  Bh  /g'  -e   's/^108 /  Hs  /g'  -e   's/^109 /  Mt  /g'

 

カテゴリー: 未分類 | コメントをどうぞ

Macportsでatlasは入れんでもいいかも。

新しいMacBook Pro(ただし最新版じゃない)を手に入れたのでLapackだのBlasだのatlasだの線形演算のライブラリをMacportsで入れなくちゃと思ってた。調べるとatlasを入れなくても、Accelerate.framworkで最適化されたBLASが提供されているのか… 知らなかった。勉強になった。

これらのリンクが参考になった。感謝。
http://decafish.blog.so-net.ne.jp/2013-03-19-1

http://www.onthelambda.com/2013/12/26/the-performance-gains-from-switching-rs-linear-algebra-libraries/

http://joker.hatenablog.com/entry/2014/11/20/013332

以下引用:

cd /Library/Frameworks/R.framework/Resources/lib
ln -sf /System/Library/Frameworks/Accelerate.framework/Versions/Current/Frameworks/vecLib.framework/Versions/Current/libBLAS.dylib libRblas.dylib

とすることで適用される。

 

 

 

 

カテゴリー: 未分類 | コメントをどうぞ

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;
}

カテゴリー: 未分類 | コメントをどうぞ

Gaussianで複合的な基底関数を使う例。

%nprocshared=16
%chk=mol.chk
%mem=24GB
#p ub3lyp/Gen scf=(maxcycle=400)
gfinput pop=full

iron metal porphyrin-bridge-cupper complex

1 1
Cu   8.273873      1.157983      -0.530262
Cu   6.679951      -0.665822      -0.293468
C   9.161808      5.035061      -2.179335
C   8.254118      4.039180      -2.080807

H   -3.910291      8.883800      2.573516
H   8.594102      1.343646      5.470718
C H N O S 0
6-31G(d,p)
****
Cu Fe 0
LANL2DZ
****

 

カテゴリー: 未分類 | コメントをどうぞ

行列をある1次元配列の順序をもとに並べ替えたい時の対処法

タイトルも何かイマイチ書きたいことを表現しきれていない気がするけど、例えば固有値問題を解いて、対角線上に固有値のある行列とそれに対応する固有ベクトルの行列を[V,Eg]=eig(H)で得ても、対角の固有値が昇順(あるいは降順)にソートされていなくって困った時に次のようにすればいいというメモ。改善の余地はあると思うけどとりあえず。

(1) 次のような関数を作る。戻り値が複数あるときは[]でくくるのかな?この場合はv1とdに帰ってくる。
% solve eigen-value problem and sort eigenvalue and vectors afterward
function [v1,d] = eigsort(x)
[v,d]   = eig (x);
[dd,ix] = sort (diag (d));
v1      = v(:,ix) ;
d=diag(sort (diag (d)));
endfunction

(2) 例えばヒュッケルとかで適当にハミルトニアンを作る。

H=[ 0 1 0 0; 1 0 1 0; 0 01 0 1; 0 0 1 0;]

(3) 先ほど(1)で定義した関数を使う。
[V, Eg] = eigsort(H);

 

追記1: 複素数の時は以下のように修正した。上はエルミート限定っぽい。”ascend”と”descend”で昇順降順を入れ替えれる。

%%%%%%%%%%%%%%
function [v1,d] = eigsort(x)
[v,d]   = eig (x);
[dd,ix] =  sort (diag(real(d)),’descend’  )
v1      = v(:, ix) ;
d=diag(sort (diag (real(d))));
endfunction

追記2:
追記1の最後の行を間違ってるかも。
d=diag(    sum(d(:,ix))’    )が正しいっぽい?

 

カテゴリー: 未分類 | コメントをどうぞ

メビウスの輪をOctaveで描いてみた。

octaveでメビウスの輪をプロットしてみた。-0.2<r<0.2, 0<t<πの2つの媒介変数表示で
x= (cos 2t)*(r cos t +2 )
y= (sin 2t)*(r cos t +2 )
z= r sint
なだけなんだけど、どうもこのドット演算子の使い方が今ひとつわからない。行き当たりばったりでとりあえずプロットできたという感じ。
>> ra=linspace(-0.2,0.2,21);
>> ta=linspace(0,pi,31);
>> [r,t]=meshgrid(ra,ta);
>> x=cos(2*t).*(r.*cos(t).+2);
>> y=sin(2*t).*(r.*cos(t).+2);
>> z=r.*sin(t)
>> surface(x,y,z)
>> axis ([ -2.5 2.5 -2.5 2.5 -1 1])

Screen Shot 2016-10-10 at 10.43.38 PM.png

カテゴリー: 未分類 | コメントをどうぞ

ImageMagickのconvertでtrimがうまくいかなかった時の対処法。

PNGファイルの余白を除去しようと-trimを使ったけど、うまくいかなかった。そして背景が透明だったのでついでに白で埋めたいと思った。次のコマンドでうまくいった。

convert -trim -flatten +repage image-xxx.png image.png

 

 

カテゴリー: 未分類 | コメントをどうぞ