ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/samples/Madelung/dipoles/buildDipolarArray
Revision: 1889
Committed: Tue Jun 18 17:54:20 2013 UTC (12 years, 2 months ago) by gezelter
File size: 6870 byte(s)
Log Message:
Working on madelung energy module for multipoles

File Contents

# Content
1 #! /usr/bin/env python
2
3 """Dipolar Lattice Builder
4
5 Creates cubic lattices of dipoles to test the dipole-dipole
6 interaction code.
7
8 Usage: buildDipolarArray
9
10 Options:
11 -h, --help show this help
12 -a, --array-type-A use array type "A" (default)
13 -b, --array-type-B use array type "B"
14 -l, --lattice=... use the specified lattice ( SC, FCC, or BCC )
15 -d, --direction=... use dipole orientation (001, 111, or 011)
16 -c, --constant=... use the specified lattice constant
17 -n use the specified number of unit cells
18 -o, --output-file=... use specified output (.xyz) file
19
20 Type "A" arrays have nearest neighbor strings of antiparallel dipoles.
21
22 Type "B" arrays have nearest neighbor strings of antiparallel dipoles
23 if the dipoles are contained in a plane perpendicular to the dipole
24 direction that passes through the dipole.
25
26 Example:
27 buildDipolarArray -a -l fcc -d 001 -c 5 -n 3 -o A_fcc_001.xyz
28
29 """
30
31 __author__ = "Dan Gezelter (gezelter@nd.edu)"
32 __version__ = "$Revision: 1639 $"
33 __date__ = "$Date: 2011-09-24 16:18:07 -0400 (Sat, 24 Sep 2011) $"
34
35 __copyright__ = "Copyright (c) 2013 by the University of Notre Dame"
36 __license__ = "OpenMD"
37
38 import sys
39 import getopt
40 import string
41 import math
42 import numpy
43
44 def usage():
45 print __doc__
46
47 def createLattice(latticeType, latticeNumber, latticeConstant, dipoleDirection, doTypeA, outputFileName):
48 a1 = numpy.array([1.0,0.0,0.0])
49 a2 = numpy.array([0.0,1.0,0.0])
50 a3 = numpy.array([0.0,0.0,1.0])
51
52 sc = numpy.array([[0.0,0.0,0.0]])
53
54 bcc = numpy.array([[0.0,0.0,0.0],
55 [0.5,0.5,0.5]])
56
57 fcc = numpy.array([[0.0,0.0,0.0],
58 [0.5,0.5,0.0],
59 [0.0,0.5,0.5],
60 [0.5,0.0,0.5]])
61
62 dipole001 = numpy.array([0.0,0.0,1.0])
63 dipole011 = numpy.array([0.0,1.0/math.sqrt(2.0),1.0/math.sqrt(2.0)])
64 dipole111 = numpy.array([1.0/math.sqrt(3.0),1.0/math.sqrt(3.0),1.0/math.sqrt(3.0)])
65
66 if (latticeType.lower() == 'sc'):
67 basis_atoms = sc
68 elif(latticeType.lower() == 'fcc'):
69 basis_atoms = fcc
70 elif(latticeType.lower() == 'bcc'):
71 basis_atoms = bcc
72 else:
73 print "unknown lattice type"
74 print __doc__
75
76 ijklvector = numpy.array([1,1,1,1])
77
78 if (dipoleDirection.lower() == '001'):
79 basis_dipoles = numpy.array([dipole001, -dipole001])
80 elif(dipoleDirection.lower() == '011'):
81 basis_dipoles = numpy.array([dipole011, -dipole011])
82 elif(dipoleDirection.lower() == '111'):
83 basis_dipoles = numpy.array([dipole111, -dipole111])
84 else:
85 print "unknown lattice type"
86 print __doc__
87
88 if (not doTypeA):
89 if (latticeType.lower() == 'sc'):
90 if (dipoleDirection.lower() == '001'):
91 ijkl = numpy.array([1,1,0,0])
92 elif(dipoleDirection.lower() == '011'):
93 ijkl = numpy.array([1,1,-1,0])
94 elif(dipoleDirection.lower() == '111'):
95 ijkl = numpy.array([1,1,-2,0])
96 elif(latticeType.lower() == 'fcc'):
97 if (dipoleDirection.lower() == '001'):
98 ijkl = numpy.array([1,1,0,1])
99 elif(dipoleDirection.lower() == '011'):
100 ijkl = numpy.array([1,1,-1,2])
101 elif(dipoleDirection.lower() == '111'):
102 ijkl = numpy.array([1,1,-2,1])
103 elif(latticeType.lower() == 'bcc'):
104 if (dipoleDirection.lower() == '001'):
105 ijkl = numpy.array([1,1,0,1])
106 elif(dipoleDirection.lower() == '011'):
107 ijkl = numpy.array([1,1,-1,1])
108 elif(dipoleDirection.lower() == '111'):
109 ijkl = numpy.array([1,1,-2,0])
110 else:
111 if (latticeType.lower() == 'sc'):
112 ijkl = numpy.array([1,1,1,0])
113 elif(latticeType.lower() == 'fcc'):
114 ijkl = numpy.array([1,1,1,1])
115 elif(latticeType.lower() == 'bcc'):
116 ijkl = numpy.array([0,0,0,1])
117
118 bravais_lattice = []
119 bravais_dipoles = []
120
121 for i in range(latticeNumber):
122 for j in range(latticeNumber):
123 for k in range(latticeNumber):
124 for l in range(len(basis_atoms)):
125 bravais_lattice.append(i*a1 + j*a2 + k*a3 + basis_atoms[l])
126 tester = numpy.array([i,j,k,l])
127 flip = numpy.dot(tester, ijkl) % 2
128 print flip
129 bravais_dipoles.append(basis_dipoles[flip])
130
131 outputFile = open(outputFileName,'w')
132
133 outputFile.write(repr(len(bravais_lattice)) + '\n')
134 Hxx = latticeConstant * latticeNumber
135 Hyy = latticeConstant * latticeNumber
136 Hzz = latticeConstant * latticeNumber
137
138 outputFile.write('Hmat: {{%d, 0, 0}, {0, %d, 0}, {0, 0, %d}}\n' % (Hxx, Hyy, Hzz))
139
140 for i in range(len(bravais_lattice)):
141 xcart = latticeConstant*(bravais_lattice[i][0])
142 ycart = latticeConstant*(bravais_lattice[i][1])
143 zcart = latticeConstant*(bravais_lattice[i][2])
144 dx = bravais_dipoles[i][0]
145 dy = bravais_dipoles[i][1]
146 dz = bravais_dipoles[i][2]
147 outputFile.write('Ar ' + repr(xcart) + ' ' + repr(ycart) + ' ' + repr(zcart) + ' ' + repr(dx) + ' ' + repr(dy) + ' ' + repr(dz) + '\n')
148
149 outputFile.close()
150
151 def main(argv):
152
153 doTypeA = True
154 haveOutputFileName = False
155 latticeType = "fcc"
156 dipoleDirection = "001"
157 latticeNumber = 3
158 latticeConstant = 4
159 try:
160 opts, args = getopt.getopt(argv, "habl:d:c:n:o:", ["help","array-type-A", "array-type-B", "lattice=" "direction=", "constant=", "output-file="])
161 except getopt.GetoptError:
162 usage()
163 sys.exit(2)
164 for opt, arg in opts:
165 if opt in ("-h", "--help"):
166 usage()
167 sys.exit()
168 elif opt in ("-a", "--array-type-A"):
169 doTypeA = True
170 elif opt in ("-b", "--array-type-B"):
171 doTypeA = False
172 elif opt in ("-l", "--lattice"):
173 latticeType = arg
174 elif opt in ("-d", "--direction"):
175 dipoleDirection = arg
176 elif opt in ("-c", "--constant"):
177 latticeConstant = float(arg)
178 elif opt in ("-n"):
179 latticeNumber = int(arg)
180 elif opt in ("-o", "--output-file"):
181 outputFileName = arg
182 haveOutputFileName = True
183 if (not haveOutputFileName):
184 usage()
185 print "No output file was specified"
186 sys.exit()
187
188 createLattice(latticeType, latticeNumber, latticeConstant, dipoleDirection, doTypeA, outputFileName);
189
190 if __name__ == "__main__":
191 if len(sys.argv) == 1:
192 usage()
193 sys.exit()
194 main(sys.argv[1:])

Properties

Name Value
svn:executable *