ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/samples/Madelung/dipoles/buildDipolarArray
Revision: 1905
Committed: Wed Jul 17 20:39:58 2013 UTC (12 years, 1 month ago) by gezelter
File size: 14193 byte(s)
Log Message:
fixing Luttinger Tisza stuff for Madan, fixing Inversions and wild Cards for James

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

Properties

Name Value
svn:executable *