# How to start with biaxial tensile strain and compressive strain in TMDCs using VASP?

Good morning!

I would like to test the biaxial tensile strain and compressive strain in a hexagonal material from the TMDC group (for example HfS2). How to get it into VASP? What are the best ways to optimize geometry? Can I fall into some traps?

HfS2
1.0
3.6389749050         0.0000000000         0.0000000000
-1.8194881106         3.1514448821         0.0000000000
0.0000000000         0.0000000000        22.8899211884
Hf    S
1    2
Direct
0.000000000         0.000000000         0.500000000
0.666666985         0.333332986         0.563098013
0.333332986         0.666666985         0.436901987


Could you help me with this as example?

• +1. Thank you Milosz for contributing your very nice question here. Could you please put the contents of your screenshot, into a code block instead of doing a screenshot? What you have done is not the right way to present your question here. Commented Sep 1, 2020 at 18:58

## 2 Answers

How to get it into VASP?

The key point is to generate different POSCAR files. Once you prepare input files, you can perform the calculation with VASP. I will assume you are considering monolayer T-phase HfS2 and show how to generate compressive and tensile structures.

• Initial structure:[HfS2.vasp]

HfS2
1.0
3.6389749050         0.0000000000         0.0000000000
-1.8194881106         3.1514448821         0.0000000000
0.0000000000         0.0000000000        22.8899211884
Hf    S
1    2
Direct
0.000000000         0.000000000         0.500000000
0.666666985         0.333332986         0.563098013
0.333332986         0.666666985         0.436901987

• Then you can apply the biaxial strain for this structure by changing the length of $$\vec{a}$$ and $$\vec{b}$$ at the same time. In detail, you can do this with the following python script:

import numpy as np

def bistrain(path1,path2,strain):

with open(path1,'r') as f1:
lines = f1.readlines()

lattice = np.zeros((3,3))
for i in range(3):lattice[i,:]=list(map(float,lines[2+i].strip().split()))

lattice[0:2,0:2]=lattice[0:2,0:2]*strain

with open(path2+'strain_'+str(strain)+'.vasp','w') as f2:
f2.write(lines[0])
f2.write(lines[1])

for j in range(3):f2.write("%20.16f  %20.16f  %20.16f" %(lattice[j,0],lattice[j,1],lattice[j,2])+'\n')

for k in range(5,len(lines)):f2.write(lines[k])

#==================================================================
path1='./HfS2.vasp'
path2='./'
#strain=0.99   ## compressive strain
strain=1.01    ## tensile strain
bistrain(path1,path2,strain)


What are the best ways to optimize geometry?

When the POSCAR is prepared, you can relax the structure referring to this answer.

Can I fall into some traps?

For the initial structure, you should use fractional coordinates rather than Cartesian coordinates.

May it hopes.

• can you please explain the script. its not generating any output Commented Nov 10, 2022 at 8:44
• It looks like this was intended to be a comment on one of the answers. When you have sufficient rep, you can add a comment. For now, I'm moving this onto Jack's answer. If you don't know Python well, you might not realize his script is a function that you need to explicitly call, rather than just running the file containing it.
– Tyberius
Commented Nov 10, 2022 at 13:05

Just to add to jack's answer. This is easily implemented in ASE.

Here is the same example which creates biaxial strain POSCARs from -5% to 5% strain in 1% increments.

from ase.io import read
from os import makedirs
import numpy as np

for strain in np.arange(0.95, 1.05, 0.01):
atoms = read("POSCAR")
atoms.cell[0:2, 0:2] *= strain
atoms.positions[:, 0:2] *= strain
makedirs("biaxial_{}".format(strain))
atoms.write("biaxial_{}/POSCAR".format(strain))