ABINIT Discussion Forums

Code: Select all

#!/usr/bin/env python

import os,string,sys,math

#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=string.split(line)[1:]
                atomic_positions.append(coord)
        if(string.find(line,"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

        #Build input file for cut3d
        f=open('sphere','w')
        f.write("tspin_2o_DEN\n")
        f.write("1\n3\n1\n1\n")
        f.write("0.0 0.0 0.0\n")
        x0=string.atof(atomic_positions[iatom][0])-cube_side/2.0
        y0=string.atof(atomic_positions[iatom][1])-cube_side/2.0
        z0=string.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(`x`+' '+`y`+' '+`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:
                        npts_integral=npts_integral+1
                        g.write(string.split(line)[4]+'\n')
                        sum=sum+string.atof(string.split(line)[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