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
|
#!/usr/bin/env python
from __future__ import print_function, division
import os,string,sys,math
from locale import atof
#CUT3D="../../../../src/98_main/cut3d"
CUT3D="cut3d"
#First find the coordinates of atoms
f=open('first_round','w')
f.write("tspin_2o_DEN\n")
f.write("1\n3\n0\n0\n")
f.close()
start_atomic=0
atomic_positions=[]
for line in os.popen(CUT3D+"< first_round").readlines():
if(start_atomic==1):
if(len(line)<3):
start_atomic=0
break
coord=line.split()[1:]
atomic_positions.append(coord)
#if(string.find(line,"Atomic positions")>0):
if(line.find("Atomic positions")>0):
start_atomic=1
print("Number of atoms = ", len(atomic_positions))
print("Atomic coordinates")
print(atomic_positions)
num_atoms=len(atomic_positions)
#Define cube
cube_side=3.3
cube_step=0.1
radius=cube_side/2.0
npts=int(cube_side/cube_step)
print("number of integration points:",npts)
#Loop on atoms
g=open('data','w')
integral=[]
for iatom in range(num_atoms):
npts_integral=0
sum=0.0
#Build input file for cut3d
f=open('sphere','w')
f.write("tspin_2o_DEN\n")
#f.write("1\n3\n1\n1\n")
f.write("3\n1\n1\n")
f.write("0.0 0.0 0.0\n")
x0=atof(atomic_positions[iatom][0])-cube_side/2.0
y0=atof(atomic_positions[iatom][1])-cube_side/2.0
z0=atof(atomic_positions[iatom][2])-cube_side/2.0
print("Treating atom #",iatom)
for i in range(npts):
for j in range(npts):
for k in range(npts):
#print i,j,k
x=x0+i*cube_step
y=y0+j*cube_step
z=z0+k*cube_step
#if((x-x0)**2+(y-y0)**2+(z-z0)**2<radius**2):
f.write('1'+"\n"+'1'+'\n'+'1'+'\n')
f.write(repr(x)+' '+repr(y)+' '+repr(z)+'\n')
f.write('0\n')
f.close()
#Execute cut3d
for line in os.popen(CUT3D+"< sphere").readlines():
#if string.find(line,"Spin difference")!=-1:
if line.find("Spin difference") != -1:
npts_integral=npts_integral+1
g.write(line.split()[4]+'\n')
sum=sum+atof(line.split()[4])
integral.append(sum/npts_integral)
#Compute magnetic moment
vol=4.*math.pi*radius**3/4.
vol=cube_side**3
for iatom in range(num_atoms):
print("For atom",iatom,"magnetic moment",integral[iatom]*vol)
|
