Skip to content

Simulation of dipole spirals in stretched PTO membranes

First is the construction of the "dipole spiral" lattice and its atomic positions.

This is particularly important in DFT calculations. Due to DFT's high dependency on the initial configuration, it is difficult to obtain a dipole spiral simply by applying external factors to the T-phase.

The code for constructing the dipole spiral achieves this by meticulously designing the atomic positions, resulting in an ideal dipole spiral structure. The specific script is as follows, where:

  • theta represents the interlayer rotation angle.

  • phase_angle is the polarization in-plane phase angle (which can be designed according to specific requirements).

  • xyangle is the angle between the two displacement chains of Pb and Ti.

It is worth noting that the dipole spiral formed in PTO can only be stabilized under in-plane tensile strain.

In finite-temperature MD simulations, the emergence of dipole spirals is disfavored by in-plane supercell sizes that are either too small (e.g., \(1 \times 1\), \(2 \times 2\), \(3 \times 3\), \(4 \times 4\)) or too large (e.g., \(30 \times 30\)). Similarly, out-of-plane supercell dimensions that are too large or too small (e.g., 1, 2, or 30 layers) also impede their formation.

In 0K DFT simulations, aided by periodic boundary conditions, dipole spirals can be stably maintained in supercells of at least \(1 \times 1 \times N\) (where \(N > 4\)).

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
import os
import numpy as np
import sys
import random
import math

ux = 3.970
uy = 3.970
uz = 3.90122575209021534000

nx = 1
ny = 1
nz = int(sys.argv[1])
natom = 5 * nx * ny * nz
nBatoms = nx * ny * nz

lx = nx * ux
ly = ny * uy
lz = nz * uz

DPbxy = 0.31427807479841097920
DPbz = 0.0464832009
DTixy = 0.20180574293849946598
DTiz = 0.0301010151

theta = (360 / nz) * 1
theta_radians = math.radians(theta) #how many pi

phase_piece = int(sys.argv[2])
phase_angle = phase_piece * np.pi/12.0

xy_angle = int(sys.argv[3])
xyangle = xy_angle * np.pi/180.0
#delta = 0.0

f2 = open('PTO_spiral.data', 'w')
f2.write("\n")
f2.write("%d atoms\n" % natom)
f2.write(" 4 atom types\n")
f2.write(" %f %f xlo xhi\n" % (0.0, lx))
f2.write(" %f %f ylo yhi\n" % (0.0, ly))
f2.write(" %f %f zlo zhi\n" % (0.0, lz))
f2.write(" 0.0000000000 0.0000000000 0.0000000000 xy xz yz\n")
f2.write("\n")
f2.write("Atoms # atomic\n")
f2.write("\n")

f3 = open('PTO_spiral.xsf', 'w')
f3.write("CRYSTAL\n")
f3.write("PRIMVEC\n")
f3.write("%15.10f %15.10f %15.10f \n" % (lx, 0.0, 0.0))
f3.write("%15.10f %15.10f %15.10f \n" % (0.0, ly, 0.0))
f3.write("%15.10f %15.10f %15.10f \n" % (0.0, 0.0, lz))
f3.write("PRIMCOORD\n")
f3.write("   %d   %d \n" % (natom, 1))

#id mol type  x y z
#coord1,molID,typeID,charge
ID = 0
#molID = 1
for k in range(0, nz):
    for j in range(0, ny):
        for i in range(0, nx):
            # Pb
            x = 0.0 + i * ux + DPbxy * math.cos(k * theta_radians + phase_angle + xyangle)
            y = 0.0 + j * uy + DPbxy * math.sin(k * theta_radians + phase_angle + xyangle)
            z = 0.0 + k * uz + DPbz
            ID += 1
            f2.write('%d %d %12.10f %12.10f %12.10f \n' % (ID, 2, x, y, z))
            f3.write('%s  %12.10f %12.10f %12.10f \n' % ("Pb", x, y, z))


for k in range(0, nz):
    for j in range(0, ny):
        for i in range(0, nx):
            # Ti
            x = ux * 0.5 + i * ux + DTixy * math.cos(k * theta_radians + phase_angle)
            y = uy * 0.5 + j * uy + DTixy * math.sin(k * theta_radians + phase_angle)
            z = ux * 0.5 + k * uz + DTiz
            ID += 1
            f2.write('%d %d %12.10f %12.10f %12.10f \n' % (ID, 3, x, y, z))
            f3.write('%s  %12.10f %12.10f %12.10f \n' % ("Ti", x, y, z))


for k in range(0, nz):
    for j in range(0, ny):
        for i in range(0, nx):
            #O1
            x = ux * 0.5 + i * ux
            y = uy * 0.5 + j * uy
            z = ux * 0.0 + k * uz
            ID += 1
            f2.write('%d %d %12.10f %12.10f %12.10f \n' % (ID, 4, x, y, z))
            f3.write('%s  %12.10f %12.10f %12.10f \n' % ("O", x, y, z))
            #O2
            x = ux * 0.5 + i * ux
            y = uy * 0.0 + j * uy
            z = ux * 0.5 + k * uz
            ID += 1
            f2.write('%d %d %12.10f %12.10f %12.10f \n' % (ID, 4, x, y, z))
            f3.write('%s  %12.10f %12.10f %12.10f \n' % ("O", x, y, z))
            #O3
            x = ux * 0.0 + i * ux
            y = uy * 0.5 + j * uy
            z = ux * 0.5 + k * uz
            ID += 1
            f2.write('%d %d %12.10f %12.10f %12.10f \n' % (ID, 4, x, y, z))
            f3.write('%s  %12.10f %12.10f %12.10f \n' % ("O", x, y, z))


f2.write("\n")
f3.write("\n")

f2.close()
f3.close()
#print "ATOMS:", ID
print("ATOMS:", ID)