## Written by Xia Qian
import numpy as np
from scipy import optimize, integrate
import sys
import getopt
def readPOSCAR(filename='POSCAR-unitcell'):
print('Open file:', filename, 'to read (as unit cell)')
f = open(filename, 'r')
# first line
cellname=f.readline().strip()
# second line
scale = float(f.readline().strip().split()[0])
## line 3-5
lat=np.zeros((3,3),dtype=np.float64)
for i in range(3):
s = f.readline().strip().split()
for j in range(3):
lat[i, j] = float(s[j]) * scale
## line 6, namelist
namelist = f.readline().strip().split()
nTypes = len(namelist)
## line 7, numlist
numlist = [0] * nTypes
s = f.readline().strip().split()
for i in range(nTypes):
numlist[i] = int(s[i])
nAtoms = sum(numlist)
## line 8, C/K or D
s = f.readline().strip()
if s[0] == 'C' or s[0] == 'c' or s[0] == 'K' or s[0] == 'k':
cart=True
else:
cart=False
## line 9-, coordinates
fraclist = np.zeros((nAtoms, 3), dtype=np.float64)
cartlist = np.zeros((nAtoms, 3), dtype=np.float64)
for i in range(nAtoms):
s = f.readline().strip().split()
for j in range(3):
fraclist[i, j] = float(s[j])
if cart is True:
cartlist[i, :] = fraclist[i, :] * scale
fraclist[i, :] = cart2frac(cartlist[i, :], lat)
else:
cartlist[i, :] = frac2cart(fraclist[i, :], lat)
f.close()
return lat, nAtoms, namelist, numlist, cartlist, fraclist
def writePOSCAR(filename, cellname, lat, nAtoms, namelist, numlist, cartlist, fraclist, cart=True, fixed=False):
print('Open file:', filename, 'to write')
if fixed == True:
sFixed="F F F"
else:
sFixed="T T T"
f = open(filename, 'w')
f.write('%s\n' % cellname)
f.write('%19.14f\n' % 1.0)
for i in range(3):
f.write('%23.16f\t%23.16f\t%23.16f\n' % (lat[i, 0],lat[i, 1],lat[i, 2]))
for i in namelist:
f.write('%4s ' % i)
f.write('\n')
for i in numlist:
f.write('%4d ' % i)
f.write('\n')
f.write('Selective dynamics\n')
if cart is True:
f.write('Cart\n')
for i in range(nAtoms):
if fixed:
f.write('%16.8f\t%16.8f\t%16.8f\tF\tF\tF\n' % (cartlist[i,0],cartlist[i,1],cartlist[i,2]))
else:
f.write('%16.8f\t%16.8f\t%16.8f\tT\tT\tT\n' % (cartlist[i,0],cartlist[i,1],cartlist[i,2]))
else:
f.write('Direct\n')
for i in range(nAtoms):
if fixed:
f.write('%16.8f\t%16.8f\t%16.8f\tF\tF\tF\n' % (fraclist[i,0],fraclist[i,1],fraclist[i,2]))
else:
f.write('%16.8f\t%16.8f\t%16.8f\tT\tT\tT\n' % (fraclist[i,0],fraclist[i,1],fraclist[i,2]))
f.close()
def cart2frac(cart, latt):
# mt=np.transpose(latt)
mt = np.linalg.inv(latt)
mt = np.transpose(mt)
frac = np.matmul(mt, cart)
return frac
def frac2cart(frac, latt):
mt = np.transpose(latt)
cart = np.matmul(mt, frac)
return cart
# function of curv length: g(x)=\sqrt(1-y'^2)
def g(x,A,w,phi):
y=np.sqrt(1+A**2*w**2*(np.cos(w*x+phi))**2)
return y
# function f(x)=\int_0^x {g(x)dx}
def func(x,s,A,w,phi):
itgl=integrate.quad(g,0,x,(A,w,phi));
y=itgl[0]-s
return y
# get parameter A of sine curve
def fa(A,L0,L,w,phi):
#itgl=quad(@g,0,L,[],[],A,w,phi)
itgl=integrate.quad(g,0,L,(A,w,phi))
#print('Here',itgl)
y=itgl[0]-L0
return y
# expand the unit cell along a,b,c
def expandCell(dim, lat_uc, nAtoms_uc, namelist_uc, numlist_uc, cartlist_uc):
nCells=dim[0]*dim[1]*dim[2]
lat=np.zeros((3,3),dtype=np.float64)
for i in range(3):
lat[i,:]=lat_uc[i,:]*dim[i]
nAtoms=nAtoms_uc*nCells
fraclist = np.zeros((nAtoms, 3), dtype=np.float64)
cartlist = np.zeros((nAtoms, 3), dtype=np.float64)
namelist=namelist_uc
numlist=numlist_uc
for i in range(len(numlist)):
numlist[i]=numlist_uc[i]*nCells
#R=np.zeros((3, 1), dtype=np.float64)
iCount=0
for at in range(nAtoms_uc):
for i in range(dim[0]):
Ra=i*lat_uc[0,:]
for j in range(dim[1]):
Rb=j*lat_uc[1,:]
for k in range(dim[2]):
Rc=k*lat_uc[2,:]
R=Ra+Rb+Rc
#print(at,i,j,k,cartlist[at,:])
cartlist[iCount,:]=cartlist_uc[at,:]+R[:]
iCount = iCount +1
fraclist[at, :] = cart2frac(cartlist[at, :], lat)
return lat, nAtoms, namelist, numlist, cartlist, fraclist
# shrink the cell along a
def shrinkCell(eta, lat, nAtoms, cartlist, fraclist):
ddL=1.0-eta/100.0
lat[0,:]=lat[0,:]*ddL # Shrink only along direction a
fraclistShrink = np.zeros((nAtoms, 3), dtype=np.float64)
cartlistShrink = np.zeros((nAtoms, 3), dtype=np.float64)
for at in range(nAtoms):
cartlistShrink[at,0]=cartlist[at,0]*ddL
cartlistShrink[at,1]=cartlist[at,1]
cartlistShrink[at,2]=cartlist[at,2]
fraclistShrink[at,:] = cart2frac(cartlist[at,:], lat)
return lat, cartlistShrink, fraclistShrink
# main function starts here:
# handle the input command line:
def main():
NC=10
centerAtom=49
print('Unit cell has been expanded by',NC,'times!')
## read POSCAR file (unit cell)
lat_uc, nAtoms_uc, namelist_uc, numlist_uc, cartlist_uc, fraclist_uc = readPOSCAR()
## expand the cell
lat, nAtoms, namelist, numlist, cartlist, fraclist = expandCell([NC,1,1], lat_uc, nAtoms_uc, namelist_uc, numlist_uc, cartlist_uc)
## write flat
cellname='flat_NC'+str(NC)+'_eta'
foFlatName='POSCAR_flat_NC'+str(NC)+'_eta'+'.vasp'
writePOSCAR(foFlatName,cellname, lat, nAtoms, namelist, numlist, cartlist, fraclist, cart=True, fixed=False)
L0=lat[0,0]
print(L0)
radius=L0/2/np.pi
print(radius)
print(cartlist[0][1])
print(nAtoms)
center=cartlist[centerAtom][2]
print(center)
for i in range(0,nAtoms):
theta=cartlist[i][0]/radius
# print(theta/2/np.pi*360)
aa=(cartlist[i][2]-center)*np.sin(theta)
bb=(cartlist[i][2]-center)*np.cos(theta)
cartlist[i][2]=radius*np.sin(theta)+ aa # z-direction
cartlist[i][0]=radius*np.cos(theta)+ bb # x-direction
cellname='tube_NC'+str(NC)+'_eta'
foFlatName='POSCAR_tube_NC'+str(NC)+'_eta'+'.vasp'
writePOSCAR(foFlatName,cellname, lat, nAtoms, namelist, numlist, cartlist, fraclist, cart=True, fixed=False)
if __name__ == "__main__":
main()
#!/public/software/anaconda3/bin/python
## flexure with two-electrodes fixed
## Files to read: POSCAR (unit cell)
## Usage: flex.py NC eta
## Here NC is the period to expand the unit cell, e.g. 20 means 20 unit cells in the central area
## and eta is compressive ratio, e.g., 5 means 5%
## Example: flex.py 20 5
import numpy as np
from scipy import optimize, integrate
import sys
import getopt
def readPOSCAR(filename='POSCAR-unitcell'):
print('Open file:', filename, 'to read (as unit cell)')
f = open(filename, 'r')
# first line
cellname=f.readline().strip()
# second line
scale = float(f.readline().strip().split()[0])
## line 3-5
lat=np.zeros((3,3),dtype=np.float)
for i in range(3):
s = f.readline().strip().split()
for j in range(3):
lat[i, j] = float(s[j]) * scale
## line 6, namelist
namelist = f.readline().strip().split()
nTypes = len(namelist)
## line 7, numlist
numlist = [0] * nTypes
s = f.readline().strip().split()
for i in range(nTypes):
numlist[i] = int(s[i])
nAtoms = sum(numlist)
## line 8, C/K or D
s = f.readline().strip()
if s[0] == 'C' or s[0] == 'c' or s[0] == 'K' or s[0] == 'k':
cart=True
else:
cart=False
## line 9-, coordinates
fraclist = np.zeros((nAtoms, 3), dtype=np.float)
cartlist = np.zeros((nAtoms, 3), dtype=np.float)
for i in range(nAtoms):
s = f.readline().strip().split()
for j in range(3):
fraclist[i, j] = float(s[j])
if cart is True:
cartlist[i, :] = fraclist[i, :] * scale
fraclist[i, :] = cart2frac(cartlist[i, :], lat)
else:
cartlist[i, :] = frac2cart(fraclist[i, :], lat)
f.close()
return lat, nAtoms, namelist, numlist, cartlist, fraclist
def writePOSCAR(filename, cellname, lat, nAtoms, namelist, numlist, cartlist, fraclist, cart=True, fixed=False):
print('Open file:', filename, 'to write')
if fixed == True:
sFixed="F F F"
else:
sFixed="T T T"
f = open(filename, 'w')
f.write('%s\n' % cellname)
f.write('%19.14f\n' % 1.0)
for i in range(3):
f.write('%23.16f\t%23.16f\t%23.16f\n' % (lat[i, 0],lat[i, 1],lat[i, 2]))
for i in namelist:
f.write('%4s ' % i)
f.write('\n')
for i in numlist:
f.write('%4d ' % i)
f.write('\n')
f.write('Selective dynamics\n')
if cart is True:
f.write('Cart\n')
for i in range(nAtoms):
if fixed:
f.write('%16.8f\t%16.8f\t%16.8f\tF\tF\tF\n' % (cartlist[i,0],cartlist[i,1],cartlist[i,2]))
else:
f.write('%16.8f\t%16.8f\t%16.8f\tT\tT\tT\n' % (cartlist[i,0],cartlist[i,1],cartlist[i,2]))
else:
f.write('Direct\n')
for i in range(nAtoms):
if fixed:
f.write('%16.8f\t%16.8f\t%16.8f\tF\tF\tF\n' % (fraclist[i,0],fraclist[i,1],fraclist[i,2]))
else:
f.write('%16.8f\t%16.8f\t%16.8f\tT\tT\tT\n' % (fraclist[i,0],fraclist[i,1],fraclist[i,2]))
f.close()
def cart2frac(cart, latt):
# mt=np.transpose(latt)
mt = np.linalg.inv(latt)
mt = np.transpose(mt)
frac = np.matmul(mt, cart)
return frac
def frac2cart(frac, latt):
mt = np.transpose(latt)
cart = np.matmul(mt, frac)
return cart
# function of curv length: g(x)=\sqrt(1-y'^2)
def g(x,A,w,phi):
y=np.sqrt(1+A**2*w**2*(np.cos(w*x+phi))**2)
return y
# function f(x)=\int_0^x {g(x)dx}
def func(x,s,A,w,phi):
itgl=integrate.quad(g,0,x,(A,w,phi));
y=itgl[0]-s
return y
# get parameter A of sine curve
def fa(A,L0,L,w,phi):
#itgl=quad(@g,0,L,[],[],A,w,phi)
itgl=integrate.quad(g,0,L,(A,w,phi))
#print('Here',itgl)
y=itgl[0]-L0
return y
# expand the unit cell along a,b,c
def expandCell(dim, lat_uc, nAtoms_uc, namelist_uc, numlist_uc, cartlist_uc):
nCells=dim[0]*dim[1]*dim[2]
lat=np.zeros((3,3),dtype=np.float)
for i in range(3):
lat[i,:]=lat_uc[i,:]*dim[i]
nAtoms=nAtoms_uc*nCells
fraclist = np.zeros((nAtoms, 3), dtype=np.float)
cartlist = np.zeros((nAtoms, 3), dtype=np.float)
namelist=namelist_uc
numlist=numlist_uc
for i in range(len(numlist)):
numlist[i]=numlist_uc[i]*nCells
#R=np.zeros((3, 1), dtype=np.float)
iCount=0
for at in range(nAtoms_uc):
for i in range(dim[0]):
Ra=i*lat_uc[0,:]
for j in range(dim[1]):
Rb=j*lat_uc[1,:]
for k in range(dim[2]):
Rc=k*lat_uc[2,:]
R=Ra+Rb+Rc
#print(at,i,j,k,cartlist[at,:])
cartlist[iCount,:]=cartlist_uc[at,:]+R[:]
iCount = iCount +1
fraclist[at, :] = cart2frac(cartlist[at, :], lat)
return lat, nAtoms, namelist, numlist, cartlist, fraclist
# shrink the cell along a
def shrinkCell(eta, lat, nAtoms, cartlist, fraclist):
ddL=1.0-eta/100.0
lat[0,:]=lat[0,:]*ddL # Shrink only along direction a
fraclistShrink = np.zeros((nAtoms, 3), dtype=np.float)
cartlistShrink = np.zeros((nAtoms, 3), dtype=np.float)
for at in range(nAtoms):
cartlistShrink[at,0]=cartlist[at,0]*ddL
cartlistShrink[at,1]=cartlist[at,1]
cartlistShrink[at,2]=cartlist[at,2]
fraclistShrink[at,:] = cart2frac(cartlist[at,:], lat)
return lat, cartlistShrink, fraclistShrink
# main function starts here:
# handle the input command line:
def main():
NC=22
eta=10
# eta=None
# try:
# opts, args = getopt.getopt(argv,"he:",["eta="])
# except getopt.GetoptError:
# print('buckling.py -e <a positive value>')
# sys.exit(2)
# for opt,arg in opts:
# if opt == '-h':
# print('buckling.py -e <a positive value>')
# sys.exit()
# elif opt in ("-e","--eta"):
# eta=arg
print('Unit cell has been expanded by',NC,'times!')
#eta=10
ddL=1.0-eta/100.0
print(ddL,'% to shrink along a')
## read POSCAR file (unit cell)
lat_uc, nAtoms_uc, namelist_uc, numlist_uc, cartlist_uc, fraclist_uc = readPOSCAR()
## expand the cell
lat, nAtoms, namelist, numlist, cartlist, fraclist = expandCell([NC,1,1], lat_uc, nAtoms_uc, namelist_uc, numlist_uc, cartlist_uc)
L0=lat[0,0]
## shrink the lattice along a
lat, cartlistFlat, fraclistFlat=shrinkCell(eta, lat, nAtoms, cartlist, fraclist)
## write flat
cellname='flat_NC'+str(NC)+'_eta'+str(eta)
foFlatName='POSCAR_flat_NC'+str(NC)+'_eta'+str(eta)+'.vasp'
writePOSCAR(foFlatName,cellname, lat, nAtoms, namelist, numlist, cartlistFlat, fraclistFlat, cart=True, fixed=False)
## get parameters of sine curve
L=L0*ddL
w=2.0*np.pi/L
phi=0.0-np.pi/2.0
### get A
A=optimize.newton(fa,0,None,(L0,L,w,phi))
A=abs(A)
## write buckled
cartlistBend=np.zeros((nAtoms, 3), dtype=np.float)
fraclistBend=np.zeros((nAtoms, 3), dtype=np.float)
for i in range(nAtoms):
s=cartlist[i,0]
x=optimize.newton(func,0,None,(s,A,w,phi))
cartlistBend[i,0]=x
cartlistBend[i,1]=cartlist[i,1]
cartlistBend[i,2]=cartlist[i,2]+A*np.sin(w*x+phi)+A
cellname='bend_NC'+str(NC)+'_eta'+str(eta)
foBendName='POSCAR_bend_NC'+str(NC)+'_eta'+str(eta)+'.vasp'
writePOSCAR(foBendName,cellname, lat, nAtoms, namelist, numlist, cartlistBend, fraclistBend, cart=True, fixed=False)
if __name__ == "__main__":
main()
No Comments
Leave a comment Cancel