Idpp method for creating initial NEB path

IDPP method is introduced by Hannes Jónsson in his paper “Improved initial guess for minimum energy path calculations”

#!/usr/bin/python
#lipai@mail.ustc.edu.cn
import numpy as np
import os
import time
step_init=0.0001
images=input("input num of images: ")
ininame=raw_input("ini structure: ")
finname=raw_input("fin structure: ")

#read ini structure
fileopen=open(ininame,'r')
ini_data=fileopen.readlines()

head=ini_data[:9]
atom_num=sum(map(int,head[6].split()))

ini_data=ini_data[9:9+atom_num]
tmp=[]
fix=[]
for i in range(atom_num):
    tmp.append(map(float,ini_data[i].split()[0:3]))
    if(ini_data[i].split()[4]=='F'): fix.append(1)
    else: fix.append(0)
pos_a=np.array(tmp)

#read fin stucture
fileopen=open(finname,'r')
fin_data=fileopen.readlines()
fin_data=fin_data[9:9+atom_num]
tmp=[]
for i in fin_data:
    tmp.append(map(float,i.split()[0:3]))
pos_b=np.array(tmp)

#correction of periodic boundary condition
for i in range(atom_num):
    for j in range(3):
            if(pos_a[i][j]-pos_b[i][j]>0.5): pos_a[i][j]-=1
            if(pos_a[i][j]-pos_b[i][j]<-0.5): pos_b[i][j]-=1

#get dist matrix and linear interpolation
dist_a=np.zeros([atom_num,atom_num])
dist_b=np.zeros([atom_num,atom_num])
for i in range(atom_num):
    for j in range(atom_num):
            tmp_a=0
            tmp_b=0
            for k in range(3):
                    tmp_a+=(pos_a[i][k]-pos_a[j][k])**2
                    tmp_b+=(pos_b[i][k]-pos_b[j][k])**2
            dist_a[i,j]=np.sqrt(tmp_a)
            dist_b[i,j]=np.sqrt(tmp_b)
dist_im=np.zeros([images,atom_num,atom_num])
pos_im=np.zeros([images,atom_num,3])
for i in range(images):
    dist_im[i]=dist_a+(i+1.0)*(dist_b-dist_a)/(images+1.0)
    pos_im[i]=pos_a+(i+1.0)*(pos_b-pos_a)/(images+1.0)

#optimization using steepest descent method
pos_tmp=np.zeros([atom_num,3])
dist_tmp=np.zeros([atom_num,atom_num])
s0=np.zeros(images)
s1=np.zeros(images)
flag=np.zeros(images)
for im in range(images):
    if(flag[im]==1): continue
    step=step_init
    print "generate image "+str(im+1)
    loop=0
    while(1):
            for i in range(atom_num):  #get dist_tmp
                    for j in range(atom_num):
                            if(i==j):
                                    dist_tmp[i,j]=10
                            else:
                                    tmp=0
                                    for k in range(3):
                                            tmp+=(pos_im[im][i][k]-pos_im[im][j][k])**2
                                    dist_tmp[i,j]=np.sqrt(tmp)
            for i in range(atom_num):
                    for sigma in range(3):
                            grad=0
                            for j in range(atom_num):
                                    if(j!=i):
                                            grad+=2*(dist_im[im][i][j]-dist_tmp[i][j])*(pos_im[im][i][sigma]-pos_im[im][j][sigma])\
                                                            *(2*dist_im[im][i][j]-dist_tmp[i][j])/dist_tmp[i,j]**5
                            pos_tmp[i][sigma]=pos_im[im][i][sigma]+step*grad
            pos_im[im]=pos_tmp
            #judge convergence
            s0[im]=s1[im]
            s1[im]=0
            for i in range(atom_num):
                    for j in range(i):
                            s1[im]+=(dist_im[im][i][j]-dist_tmp[i][j])**2/dist_tmp[i][j]**4
            loop+=1
            print "loop:%d"%loop
            if(abs(s0[im]-s1[im])<0.01):
                    print "Converged!"
                    flag[im]=1
                    break
            if(loop>1 and s1[im]>s0[im]): step=step/3

#mkdir and generate poscar files
if(images+1<10): num='0'+str(images+1)
else:       num=str(images+1)
os.system("mkdir 00")
os.system("cp p1 00/POSCAR")
os.system("mkdir "+num)
os.system("cp p2 "+num+"/POSCAR")
for i in range(images):
    if(i+1<10): num='0'+str(i+1)
    else:       num=str(i+1)
    os.system("mkdir "+num)
    data=pos_im[i].tolist()
    filename=num+"/POSCAR"
    f=file(filename,"a+")
    f.writelines(head)
    for j in range(atom_num):
            line=map(str,data[j])
            line="  ".join(line)
            if(fix[j]==1): line=line+'    F  F  F\n'
            else:          line=line+'    T  T  T\n'
            f.write(line)
    f.close()