ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/samples/Madelung/dipoles/buildDipolarArray
Revision: 1898
Committed: Mon Jul 8 21:09:21 2013 UTC (12 years, 1 month ago) by gezelter
File size: 13139 byte(s)
Log Message: 
More work on the Luttinger & Tisza structure generator

File Contents

# User Rev Content
1 gezelter 1889 #! /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 gezelter 1898
47 gezelter 1889
48 gezelter 1898
49 gezelter 1889 def createLattice(latticeType, latticeNumber, latticeConstant, dipoleDirection, doTypeA, outputFileName):
50 gezelter 1898 # The following section creates 24 basic arrays from Luttinger and
51     # Tisza:
52 gezelter 1889
53 gezelter 1898 # The six unit vectors are: 3 spatial and 3 to describe the
54     # orientation of the dipole.
55    
56     e1 = numpy.array([1.0,0.0,0.0,0.0,0.0,0.0])
57     e2 = numpy.array([0.0,1.0,0.0,0.0,0.0,0.0])
58     e3 = numpy.array([0.0,0.0,1.0,0.0,0.0,0.0])
59     e4 = numpy.array([0.0,0.0,0.0,1.0,0.0,0.0])
60     e5 = numpy.array([0.0,0.0,0.0,0.0,1.0,0.0])
61     e6 = numpy.array([0.0,0.0,0.0,0.0,0.0,1.0])
62    
63     # Parameters describing the 8 basic arrays:
64     cell = numpy.array([[0,0,0],[0,0,1],[1,0,0],[0,1,0],
65     [1,1,0],[0,1,1],[1,0,1],[1,1,1]])
66    
67     X = numpy.zeros(384).reshape((8,8,6))
68     Y = numpy.zeros(384).reshape((8,8,6))
69     Z = numpy.zeros(384).reshape((8,8,6))
70     # mX, mY, and mZ arrays have dipole direction flipped
71     mX = numpy.zeros(384).reshape((8,8,6))
72     mY = numpy.zeros(384).reshape((8,8,6))
73     mZ = numpy.zeros(384).reshape((8,8,6))
74    
75     # create the 24 basic arrays using Eq. 12 in Luttinger & Tisza:
76     for i in range(8):
77     which = 0
78     for l1 in range(2):
79     for l2 in range(2):
80     for l3 in range(2):
81     lvals = numpy.array([l1,l2,l3])
82     value = math.pow(-1, numpy.dot(cell[i], lvals))
83     Xvec = (l1*e1 + l2*e2 + l3*e3) + value * e4
84     Yvec = (l1*e1 + l2*e2 + l3*e3) + value * e5
85     Zvec = (l1*e1 + l2*e2 + l3*e3) + value * e6
86     mXvec = (l1*e1 + l2*e2 + l3*e3) - value * e4
87     mYvec = (l1*e1 + l2*e2 + l3*e3) - value * e5
88     mZvec = (l1*e1 + l2*e2 + l3*e3) - value * e6
89     X[i][which] = Xvec
90     Y[i][which] = Yvec
91     Z[i][which] = Zvec
92     mX[i][which] = mXvec
93     mY[i][which] = mYvec
94     mZ[i][which] = mZvec
95     which = which + 1
96    
97     # The simple cubic array has only one site at the lattice point:
98     lp = numpy.array([0.0,0.0,0.0,0.0,0.0,0.0])
99    
100     # The body-centered cubic array also has a body-centerered site:
101     bc = numpy.array([0.5,0.5,0.5,0.0,0.0,0.0])
102    
103     # The face-centered cubic array also has 3 face-centered sites:
104     xy = numpy.array([0.5,0.5,0.0,0.0,0.0,0.0])
105     yz = numpy.array([0.0,0.5,0.5,0.0,0.0,0.0])
106     xz = numpy.array([0.5,0.0,0.5,0.0,0.0,0.0])
107    
108 gezelter 1889 sc = numpy.array([[0.0,0.0,0.0]])
109    
110     bcc = numpy.array([[0.0,0.0,0.0],
111     [0.5,0.5,0.5]])
112    
113     fcc = numpy.array([[0.0,0.0,0.0],
114     [0.5,0.5,0.0],
115     [0.0,0.5,0.5],
116     [0.5,0.0,0.5]])
117    
118 gezelter 1898 known_case = False
119     if (doTypeA):
120     if (latticeType.lower() == 'sc'):
121     basic_array = Z[5]+lp
122     known_case = True
123     if (latticeType.lower() == 'bcc'):
124     if (dipoleDirection.lower() == '001'):
125     basic_array = numpy.append(Z[1]+lp, mZ[1]+bc, axis=0)
126     known_case = True
127     if (dipoleDirection.lower() == '111'):
128     basic_array = numpy.append(Z[5]+X[7]+Y[6]+lp,
129     Z[5]+X[7]+Y[6]+bc, axis=0)
130     known_case = True
131     if (latticeType.lower() == 'fcc'):
132     if (dipoleDirection.lower() == '001'):
133     basic_array = numpy.append(Z[1]+lp, Z[1]+xy, axis=0)
134     basic_array = numpy.append(basic_array, mZ[1]+yz, axis=0)
135     basic_array = numpy.append(basic_array, mZ[1]+xz, axis=0)
136     known_case = True
137     if (dipoleDirection.lower() == '011'):
138     basic_array = numpy.append(Z[1]+Y[1]+lp, Z[1]+Y[1]+yz, axis=0)
139     basic_array = numpy.append(basic_array, mZ[1]+mY[1]+xy, axis=0)
140     basic_array = numpy.append(basic_array, mZ[1]+mY[1]+xz, axis=0)
141     known_case = True
142 gezelter 1889 else:
143     if (latticeType.lower() == 'sc'):
144 gezelter 1898 basic_array = Z[5]+lp
145     known_case = True
146     if (latticeType.lower() == 'bcc'):
147 gezelter 1889 if (dipoleDirection.lower() == '001'):
148 gezelter 1898 basic_array = numpy.append(Z[5]+lp, mZ[5]+bc, axis=0)
149     known_case = True
150     if (dipoleDirection.lower() == '111'):
151     basic_array = numpy.append(Z[5]+X[7]+Y[6]+lp,
152     Z[5]+X[7]+Y[6]+bc, axis=0)
153     known_case = True
154     if (latticeType.lower() == 'fcc'):
155 gezelter 1889 if (dipoleDirection.lower() == '001'):
156 gezelter 1898 basic_array = numpy.append(Z[1]+lp, Z[1]+yz, axis=0)
157     basic_array = numpy.append(basic_array, mZ[1]+xy, axis=0)
158     basic_array = numpy.append(basic_array, mZ[1]+xz, axis=0)
159     known_case = True
160     if (dipoleDirection.lower() == '011'):
161     basic_array = numpy.append(Z[8]+Y[8]+lp, Z[8]+Y[8]+xy, axis=0)
162     basic_array = numpy.append(basic_array, Z[8]+Y[8]+yz, axis=0)
163     basic_array = numpy.append(basic_array, mZ[8]+mY[8]+xz, axis=0)
164     known_case = True
165 gezelter 1889
166 gezelter 1898 if (not known_case):
167     print "unhandled combination of lattice and dipole direction"
168     print __doc__
169    
170     print basic_array
171    
172 gezelter 1889 bravais_lattice = []
173     for i in range(latticeNumber):
174     for j in range(latticeNumber):
175     for k in range(latticeNumber):
176 gezelter 1898 for l in range(len(basic_array)):
177     bravais_lattice.append(2*i*e1 + 2*j*e2 + 2*k*e3 + basic_array[l])
178 gezelter 1889
179 gezelter 1898 outputFile = open(outputFileName, 'w')
180    
181     outputFile.write('<OpenMD version=2>\n')
182     outputFile.write(' <MetaData>\n')
183     outputFile.write(' molecule{\n')
184     outputFile.write(' name = \"D\";\n')
185     outputFile.write(' \n')
186     outputFile.write(' atom[0]{\n')
187     outputFile.write(' type = \"D\";\n')
188     outputFile.write(' position(0.0, 0.0, 0.0);\n')
189     outputFile.write(' orientation(0.0, 0.0, 0.0);\n')
190     outputFile.write(' }\n')
191     outputFile.write(' }\n')
192     outputFile.write(' component{\n')
193     outputFile.write(' type = \"D\";\n')
194     outputFile.write(' nMol = '+ repr(len(bravais_lattice)) + ';\n')
195     outputFile.write(' }\n')
196 gezelter 1889
197 gezelter 1898 outputFile.write(' ensemble = NVE;\n')
198     outputFile.write(' forceField = \"Multipole\";\n')
199    
200     outputFile.write(' cutoffMethod = \"shifted_force\";\n')
201     outputFile.write(' electrostaticScreeningMethod = \"damped\";\n')
202    
203     outputFile.write(' cutoffRadius = 9.0;\n')
204     outputFile.write(' dampingAlpha = 0.18;\n')
205     outputFile.write(' statFileFormat = \"TIME|TOTAL_ENERGY|POTENTIAL_ENERGY|KINETIC_ENERGY|TEMPERATURE|PRESSURE|VOLUME|CONSERVED_QUANTITY|ELECTROSTATIC_POTENTIAL\";\n')
206     outputFile.write(' dt = 1.0;\n')
207     outputFile.write(' runTime = 1.0;\n')
208     outputFile.write(' sampleTime = 1.0;\n')
209     outputFile.write(' statusTime = 1.0;\n')
210     outputFile.write(' </MetaData>\n')
211     outputFile.write(' <Snapshot>\n')
212     outputFile.write(' <FrameData>\n');
213     outputFile.write(" Time: %.10g\n" % (0.0))
214 gezelter 1889
215 gezelter 1898 Hxx = 2.0 * latticeConstant * latticeNumber
216     Hyy = 2.0 * latticeConstant * latticeNumber
217     Hzz = 2.0 * latticeConstant * latticeNumber
218    
219     outputFile.write(' Hmat: {{%d, 0, 0}, {0, %d, 0}, {0, 0, %d}}\n' % (Hxx, Hyy, Hzz))
220     outputFile.write(' </FrameData>\n')
221     outputFile.write(' <StuntDoubles>\n')
222     sdFormat = 'pvqj'
223     index = 0
224    
225 gezelter 1889 for i in range(len(bravais_lattice)):
226     xcart = latticeConstant*(bravais_lattice[i][0])
227     ycart = latticeConstant*(bravais_lattice[i][1])
228     zcart = latticeConstant*(bravais_lattice[i][2])
229 gezelter 1898 dx = bravais_lattice[i][3]
230     dy = bravais_lattice[i][4]
231     dz = bravais_lattice[i][5]
232    
233     uz = numpy.array([dx, dy, dz])
234     uz = uz/numpy.linalg.norm(uz)
235    
236     uy = numpy.array([0.0, 1.0, 0.0])
237     uy = uy - uz * numpy.vdot(uy, uz) / numpy.vdot(uz, uz)
238     uy = uy/numpy.linalg.norm(uy)
239    
240     ux = numpy.cross(uy, uz)
241    
242     RotMat = [ux, uy, uz]
243    
244     q = [0.0, 0.0, 0.0, 0.0]
245    
246     # RotMat to Quat code is out of OpenMD's SquareMatrix3.hpp code:
247    
248     t = RotMat[0][0] + RotMat[1][1] + RotMat[2][2] + 1.0
249    
250     if( t > 1e-6 ):
251     s = 0.5 / math.sqrt( t )
252     q[0] = 0.25 / s
253     q[1] = (RotMat[1][2] - RotMat[2][1]) * s
254     q[2] = (RotMat[2][0] - RotMat[0][2]) * s
255     q[3] = (RotMat[0][1] - RotMat[1][0]) * s
256     else:
257     ad1 = RotMat[0][0]
258     ad2 = RotMat[1][1]
259     ad3 = RotMat[2][2]
260    
261     if( ad1 >= ad2 and ad1 >= ad3 ):
262     s = 0.5 / math.sqrt( 1.0 + RotMat[0][0] - RotMat[1][1] - RotMat[2][2] )
263     q[0] = (RotMat[1][2] - RotMat[2][1]) * s
264     q[1] = 0.25 / s
265     q[2] = (RotMat[0][1] + RotMat[1][0]) * s
266     q[3] = (RotMat[0][2] + RotMat[2][0]) * s
267     elif ( ad2 >= ad1 and ad2 >= ad3 ):
268     s = 0.5 / math.sqrt( 1.0 + RotMat[1][1] - RotMat[0][0] - RotMat[2][2] )
269     q[0] = (RotMat[2][0] - RotMat[0][2] ) * s
270     q[1] = (RotMat[0][1] + RotMat[1][0]) * s
271     q[2] = 0.25 / s
272     q[3] = (RotMat[1][2] + RotMat[2][1]) * s
273     else:
274     s = 0.5 / math.sqrt( 1.0 + RotMat[2][2] - RotMat[0][0] - RotMat[1][1] )
275     q[0] = (RotMat[0][1] - RotMat[1][0]) * s
276     q[1] = (RotMat[0][2] + RotMat[2][0]) * s
277     q[2] = (RotMat[1][2] + RotMat[2][1]) * s
278     q[3] = 0.25 / s
279    
280    
281     outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (index, sdFormat, xcart, ycart, zcart, 0.0, 0.0, 0.0, q[0], q[1], q[2], q[3], 0.0, 0.0, 0.0))
282     index = index+1
283    
284     outputFile.write(" </StuntDoubles>\n")
285     outputFile.write(" </Snapshot>\n")
286     outputFile.write("</OpenMD>\n")
287 gezelter 1889 outputFile.close()
288    
289 gezelter 1898 outputFile.close()
290    
291 gezelter 1889 def main(argv):
292    
293     doTypeA = True
294     haveOutputFileName = False
295     latticeType = "fcc"
296     dipoleDirection = "001"
297     latticeNumber = 3
298     latticeConstant = 4
299     try:
300     opts, args = getopt.getopt(argv, "habl:d:c:n:o:", ["help","array-type-A", "array-type-B", "lattice=" "direction=", "constant=", "output-file="])
301     except getopt.GetoptError:
302     usage()
303     sys.exit(2)
304     for opt, arg in opts:
305     if opt in ("-h", "--help"):
306     usage()
307     sys.exit()
308     elif opt in ("-a", "--array-type-A"):
309     doTypeA = True
310     elif opt in ("-b", "--array-type-B"):
311     doTypeA = False
312     elif opt in ("-l", "--lattice"):
313     latticeType = arg
314     elif opt in ("-d", "--direction"):
315     dipoleDirection = arg
316     elif opt in ("-c", "--constant"):
317     latticeConstant = float(arg)
318     elif opt in ("-n"):
319     latticeNumber = int(arg)
320     elif opt in ("-o", "--output-file"):
321     outputFileName = arg
322     haveOutputFileName = True
323     if (not haveOutputFileName):
324     usage()
325     print "No output file was specified"
326     sys.exit()
327    
328     createLattice(latticeType, latticeNumber, latticeConstant, dipoleDirection, doTypeA, outputFileName);
329    
330     if __name__ == "__main__":
331     if len(sys.argv) == 1:
332     usage()
333     sys.exit()
334     main(sys.argv[1:])

Properties

Name Value
svn:executable *