Tuesday, September 8, 2015

Substitution inside files by one line bash commands

Recently I searched at the internet for substituing line(s) in many files with increasing number. Here is the example how to substitute whole 3rd line by 'whatever' in the files named  1D*.in.

sed -i '3s/.*/whatever/' > ls 1D*.in 


How to add prefix to every file in the current directory: 

\( for filename in *; do mv "$filename" "'prefix_'$filename"; done; \)

Or eventually rename it:

rename 's/prefix_/a_/' *

or renaming nicer with command mv

\(for file in prefix_* do mv -i "${file}" "${file/prefix_/a_}" done \)

where option -i means 'in-situ'.

Now, let say we want to substitute a string 'cis' to 'fis' in every file
name by bash:

sed -i 's/cis/fis/'

Sed  is a stream editor.  A stream editor is  used to perform basic text transformations on an input stream (a file or input from a pipeline).  By default, sed sends its results to the screen. You can cancel this default printing with option -n.
Syntax: sed [options] 'commands' [file-to-edit] For commands use ampersands, like e.g.
sed '3d' - to delete the 3rd line
sed -n '2p' -to print the 2nd line and shut down the default printing on the screen by option -n. Otherwise it will print a line two times.
sed -n '/root/p' -to print only the matching line and shut down the default printing on the screen by option -n

sed -i --in-place[=SUFFIX]
 edit files in place (makes backup if SUFFIX supplied) - the files will be rewritten.

sed 's/to_be_replaced/replacement/' - if match, replace a string by a string with replacement string.

Regular characters for sed:

$    the end of line
^    the start of line
.*   everything on line
\n   new line break
$!   last line

To join more commands you can use either -e option
sed -i -e '3s/to_be_replaced/replacement/' -e '2d' *.m
or with semicolon
sed -i '3s/to_be_replaced/replacement/' ; '2d' *.m

My special thanks belongs to guys on Stackoverflow and creators of man pages for Linux.

Wednesday, September 2, 2015

Eigenvalue problem, 1D B-splines basis function

In quantum mechanics, to describe the particle or particles in space we use wave functions, which we can describe as the superposition of the basis functions.

Example of quantum state of the order o=5 in one point described by radius vector \(\vec{r}\) from the nucleus:

 \(\psi(r) = c_{1}\Phi_{1}(r) +c_{2} \Phi_{2}(r) + c_{3} \Phi_{3}(r) +c_{4}\Phi_{4}(r)+ c_{5}\Phi_{5}(r)                                                   \quad \quad \quad (1)\).

Weights here are \(c_{1} - c_{5}\).  If basis function are orthonormal, we get these weights in the range <0,1> as the eigenvector by solving eigenvalue problem \(\hat{H} \Psi(r)  =  E \Psi(r) \).
Basis functions describe whole space interval, one possibility to make atom/molecule basis functions in 1-D  are B-splines (Basis - splines). B-splines are defined in the restricted space compare to Slater or Gaussians orbitals.

So first, let's define the box [a b]:
a=0;
b=500;

Number of discrete steps defines vector of breakpoints:
dX=(b-a)/(N_b-1);
x_break=a:dX:b;

Number of B-splines, Nb, must be at least twice the length of the box b.
Nb=1000

B-splines are spanned on the whole box. Except of breaking points we define also knot-points. Each B-spline is starting on the each knot-points. Number of knot points = Number of B-splines. You can set up knot sequence as you like. 

In our case, knot points are equal to the breaking points except on the boundaries, where the multiplicity of knot-points is the order of B-splines  o+1. Also the degree of polynomials is defined as o+1.
In our case, the order is o=5. This order also define the number of control points for each splines between two knot-points.

Plot(1)


On the plot above Plot(1) are B-splines of the order o=5. You can see number of control points between each adjacent knot points is equal to the order of B-splines. Notice that at the boundary you have multiplicity of knot-points equals to of the o+1, therefore o+1 B-splines are raising from the boundary point. Notice that B-spline spreads among knot points which number equals to the order of B-splines.

Plot(2) 
Spreading of B-splines of order o=5 in the interval within 5 knot-points. Line  '--' highlights each knot point.

Plot(3)

Spreading of B-splines of order o=1 in the interval within 2 knot-points.  Line  '--' highlights
each knot point.

B-splines are given by recurrence relation

\[B_{i}^{o+1}(x)=\frac{x - t_{i}}{t_{i+o}-t_{i}}B_{i}^{o}(x) + (\frac{t_{i+o+1}-x}{t_{i+o+1} - t_{i+1}}) B_{i+1}^{o}(x)      \quad \quad \quad (2)\].
where \(B_{i}^{o}\) is the B-spline starting on the i-th knot-point t, of the order o.

Each B-spline has the non-zero amplitude in the interval between neighboring knot points:

\[B_{i}^{o}(x)=\bigg\{
\frac{>0}{0} \frac{\textbf{if} \quad  t_{i}\le x < t_{i+o} } {elsewise}\]

Be awared that you use only the B-splines of the highest order from the above reccurence equation to use it as the first approximation of eigenstates or other application. The order of splines must be at least o=4.

Each control point of the B-splines has its weight (see amplitude of the B-splines).
We get all control points (integration points) and its weight by the Golub-Welsch algorithm for first:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Gauss-Legendre integration points and  its weights%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%o =5, is order of each B-spline

i=1:o;
B=1./2./sqrt(1-(2.*i).^(-2));

J(1:o,1:o)=0;
for j=1:n-1
J(j,j+1)=B(j);
J(j+1,j)=B(j);
end
[V,D] = eig(J);

for j=1:o
WW(j)=2.*V(1,j).^2;

%%%And control points are given by diagonalization:
xx(j)=D(j,j);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Now we spread our controls points  with corresponding weights on the whole interval defined by the breaking points x_breaks. Please notice that each spline has 5 control points between each  two break points. 
%Nb is number of B-splines on the whole interval

for i=2:Nb
x((i-2)*o+1:(i-1)*o)=(x_break(i)-x_break(i-1))./2.*xx+(x_break(i)+x_break(i-1))./2;
W((i-2)*o+1:(i-1)*o)=WW;
end  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Some application:

Control points define the integration points (as well for differentiate interval).

Kinetic Energy Matrix Operator with our basis function B:

       Integral 1=-1./2.*dX./2.*sum(W'.*B(:,i+1,o+1).*diff_B(:,j+1,o+1));

        H(i,j)=H(i,j)+int1;

To remind you: i=j. The size of Hamiltonian matrix is Nb x Nb of the order o. Differentiate of B-splines, diff_B, is given by recurrence relation:

\[\frac{d}{dx}B_{i}^{o+1}(x) = \frac{o (o-1)}{t(i+o)-t(o)}.\frac{B_{i}^{o-1}(x)}{t(i+o-1)-t(i)} - \frac{B_{i+1}^{o-1}(x)}{t(i+o)-t(i+1)} \\
- \frac{o (o-1)}{t(i+o-1)-t(i)}(\frac{B_{i+1}^{o-1}(x)}{t(i+o)-t(i+1)} - \frac{B_{i+2}^{o-1}(x)}{t(i+o+1)-t(i+2)}) \quad \quad \quad (3)\].

Because B- splines are not orthogonal function, we need to compute the overlap matrix:

%S is overlap Matrix
          int2=dX./2.*sum(W'.*B(:,i+1,o+1).*B(:,j+1,o+1));
          S(i,j)=S(i,j)+int2;

After diagonalization, in matlab: [V,E] = eig(H,S), where V is the eigenstates vector, so we have new weights given by the vector V instead of initial W gained by Golub-Welsch algorithm for each B-spline in the box [a b].  E are the eigenvalues on diagonal.

We construct the state functions psi(r,i) as in the equation (1):

 psi(:,i)=psi(:,i)+V(j,i).*B(:,j+1,n+1)

j are the number of  the B-splines Nb. For each state  i you need Nb basis functions \[\Phi(r)]\, only some of them are non-zero. The number of the states i is given by the principal and angular quantum number (in this 1D case magnetic quantum number m=0). The size of V(j,j) is Nb x Nb (as large as the Hamiltonian), maximum number of the states is Nb spread on Nb points.

To summarize basic charecteristics of B-splines as the state function:

Advantages:
  • sparse Hamiltonian matrix with (Nb-1)*(2*o+1)= non-zero elements. Very compact Hamiltonian is the excellent preconditions to use it in the time propagator to solve Time Dependent Schroedinger Equation (TDSE).         
  • low RAM consumption
  • easily diagonalisable


 Disdvantages:
  • Not orthogonal, need to compute overlap matrix with make (2*(o+1)*Nb) nonzero elements
  • They are not eigenfunctions of Hamiltonian. However wavefunctions at the small distance r from nucleus are the same as the analytical solution for hydrogen, which is showed at the plot below.

Plot(4) -Energy spectrum of the Hamiltonian in B-spline basis.


Plot(5)





     For s-states green color radial wavefunctions are analytical sulutions for hydrogen (Slater-type function from
http://en.citizendium.org/wiki/Hydrogen-like_atom),  blue color wavefunctions are generated in the basis of B-splines.







For p-states red color radial wavefunctions are analytical solutions for hydrogen (Slater-type function from http://en.citizendium.org/wiki/Hydrogen-like_atom), black color wavefunctions are generated in the basis of B-splines.

To make 3D wave function $\Psi(r,\Theta,\Phi)$, multiply $\Psi(r)$ with Legendre polynomial.

Wednesday, August 12, 2015

My favourite editor - Vim


Default file: /usr/share/vim/vimfiles/archlinux.vim - setting only a few sane defaults, this file is upragading every time you upgrade/change version of vim
Global setting for root: /etc/vimrc
Adjusted for particular user: in home directory you can find .vimrc to config vim. Here are some commands for vim to manage fortran files.



" To get rid off annoying space/tab background fill.
 let fortran_have_tabs=1

" To have free source form, the same syntax in .f, .f90,  .f95, use it when you use -ffree-form as the tag in Fortran. This should be prior to the syntax comman
 let fortran_free_source=1

" Do loops will be intended
 let fortran_do_enddo=1

" The syntax coloring will be more precise but slower.  In particular,
" statement labels used in do, goto and arithmetic if statements will be
" recognized, as will construct names at the end of a do, if, select or forall construct
let fortran_more_precise=1

" Enable syntax highlighting
syntax enable

" Colorscheme
colorscheme desert
set background=dark

Explanation and more commands you can find in vim :help fortran