[ViewVC] Annotation of: OpenMD/trunk/samples/Madelung/dipoles/buildDipolarArray

#! /usr/bin/env python

"""Dipolar Lattice Builder

Creates cubic lattices of dipoles to test the dipole-dipole
interaction code.  

Usage: buildDipolarArray 

Options:
  -h, --help              show this help
  -x, --array-type-X      use one of the basic "X" arrays
  -y, --array-type-Y      use one of the basic "Y" arrays
  -z, --array-type-Z      use one of the basic "Z" arrays
  -a, --array-type-A      use array type "A" (default)
  -b, --array-type-B      use array type "B"
  -l, --lattice=...       use the specified lattice ( SC, FCC, or BCC )
  -d, --direction=...     use dipole orientation (001, 111, or 011)
  -c, --constant=...      use the specified lattice constant
  -n                      use the specified number of unit cells
  -o, --output-file=...   use specified output (.xyz) file

Type "A" arrays have nearest neighbor strings of antiparallel dipoles.

Type "B" arrays have nearest neighbor strings of antiparallel dipoles
if the dipoles are contained in a plane perpendicular to the dipole
direction that passes through the dipole.

Example:
   buildDipolarArray -a -l fcc -d 001 -c 5 -n 3 -o A_fcc_001.xyz

"""

__author__ = "Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Revision: 1639 $"
__date__ = "$Date: 2011-09-24 16:18:07 -0400 (Sat, 24 Sep 2011) $"

__copyright__ = "Copyright (c) 2013 by the University of Notre Dame"
__license__ = "OpenMD"

import sys
import getopt
import string
import math
import numpy

def usage():
    print __doc__

    
    
def createLattice(latticeType, latticeNumber, latticeConstant, dipoleDirection, arrayType, outputFileName):
    # The following section creates 24 basic arrays from Luttinger and
    # Tisza:

    # The six unit vectors are: 3 spatial and 3 to describe the
    # orientation of the dipole.
    
    e1 = numpy.array([1.0,0.0,0.0,0.0,0.0,0.0])
    e2 = numpy.array([0.0,1.0,0.0,0.0,0.0,0.0])
    e3 = numpy.array([0.0,0.0,1.0,0.0,0.0,0.0])
    e4 = numpy.array([0.0,0.0,0.0,1.0,0.0,0.0])
    e5 = numpy.array([0.0,0.0,0.0,0.0,1.0,0.0])
    e6 = numpy.array([0.0,0.0,0.0,0.0,0.0,1.0])

    # Parameters describing the 8 basic arrays:
    cell = numpy.array([[0,0,0],[0,0,1],[1,0,0],[0,1,0],
                        [1,1,0],[0,1,1],[1,0,1],[1,1,1]])

    X = numpy.zeros(384).reshape((8,8,6))
    Y = numpy.zeros(384).reshape((8,8,6))
    Z = numpy.zeros(384).reshape((8,8,6))
    # mX, mY, and mZ arrays have dipole direction flipped
    mX = numpy.zeros(384).reshape((8,8,6))
    mY = numpy.zeros(384).reshape((8,8,6))
    mZ = numpy.zeros(384).reshape((8,8,6))

    # create the 24 basic arrays using Eq. 12 in Luttinger & Tisza:
    for i in range(8):
        which = 0
        for l1 in range(2):
            for l2 in range(2):
                for l3 in range(2):
                    lvals = numpy.array([l1,l2,l3])
                    value = math.pow(-1, numpy.dot(cell[i], lvals))
                    Xvec = (l1*e1 + l2*e2 + l3*e3) + value * e4             
                    Yvec = (l1*e1 + l2*e2 + l3*e3) + value * e5             
                    Zvec = (l1*e1 + l2*e2 + l3*e3) + value * e6             
                    mXvec = (l1*e1 + l2*e2 + l3*e3) - value * e4             
                    mYvec = (l1*e1 + l2*e2 + l3*e3) - value * e5             
                    mZvec = (l1*e1 + l2*e2 + l3*e3) - value * e6             
                    X[i][which] = Xvec
                    Y[i][which] = Yvec
                    Z[i][which] = Zvec
                    mX[i][which] = mXvec
                    mY[i][which] = mYvec
                    mZ[i][which] = mZvec
                    which = which + 1

    # The simple cubic array has only one site at the lattice point:
    lp = numpy.array([0.0,0.0,0.0,0.0,0.0,0.0])

    # The body-centered cubic array also has a body-centerered site:
    bc = numpy.array([0.5,0.5,0.5,0.0,0.0,0.0])

    # The face-centered cubic array also has 3 face-centered sites:
    xy = numpy.array([0.5,0.5,0.0,0.0,0.0,0.0])
    yz = numpy.array([0.0,0.5,0.5,0.0,0.0,0.0])
    xz = numpy.array([0.5,0.0,0.5,0.0,0.0,0.0])

    sc =      numpy.array([[0.0,0.0,0.0]])

    bcc =     numpy.array([[0.0,0.0,0.0],
                           [0.5,0.5,0.5]])

    fcc =     numpy.array([[0.0,0.0,0.0], 
                           [0.5,0.5,0.0],
                           [0.0,0.5,0.5],
                           [0.5,0.0,0.5]])

    known_case = False

    if (arrayType == 'X'):
        if (int(latticeType)):
            which = int(latticeType) - 1
            basic_array = X[which]
            known_case = True
    if (arrayType == 'Y'):
        if (int(latticeType)):
            which = int(latticeType) - 1
            basic_array = Y[which]
            known_case = True
    if (arrayType == 'Z'):
        if (int(latticeType)):
            which = int(latticeType) - 1 
            basic_array = Z[which]
            known_case = True
    if (arrayType == 'A'):
        if (latticeType.lower() == 'sc'):
            basic_array = Z[5]+lp
            known_case = True
        if (latticeType.lower() == 'bcc'):
            if (dipoleDirection.lower() == '001'):
                basic_array = numpy.append(Z[1]+lp, mZ[1]+bc, axis=0)
                known_case = True
            if (dipoleDirection.lower() == '111'):
                basic_array = numpy.append(Z[5]+X[7]+Y[6]+lp,
                                           Z[5]+X[7]+Y[6]+bc, axis=0)
                known_case = True
        if (latticeType.lower() == 'fcc'):
            if (dipoleDirection.lower() == '001'):
                basic_array = numpy.append(Z[1]+lp,  Z[1]+xy, axis=0)
                basic_array = numpy.append(basic_array, mZ[1]+yz, axis=0) 
                basic_array = numpy.append(basic_array, mZ[1]+xz, axis=0) 
                known_case = True
            if (dipoleDirection.lower() == '011'):
                basic_array = numpy.append(Z[1]+Y[1]+lp, Z[1]+Y[1]+yz, axis=0)
                basic_array = numpy.append(basic_array, mZ[1]+mY[1]+xy, axis=0)
                basic_array = numpy.append(basic_array, mZ[1]+mY[1]+xz, axis=0) 
                known_case = True
    else:
        if (latticeType.lower() == 'sc'):
            basic_array = Z[5]+lp
            known_case = True
        if (latticeType.lower() == 'bcc'):
            if (dipoleDirection.lower() == '001'):
                basic_array = numpy.append(Z[5]+lp, mZ[5]+bc, axis=0)
                known_case = True
            if (dipoleDirection.lower() == '111'):
                basic_array = numpy.append(Z[5]+X[7]+Y[6]+lp,
                                           Z[5]+X[7]+Y[6]+bc, axis=0)
                known_case = True
        if (latticeType.lower() == 'fcc'):
            if (dipoleDirection.lower() == '001'):
                basic_array = numpy.append(Z[1]+lp,  Z[1]+yz, axis=0)
                basic_array = numpy.append(basic_array, mZ[1]+xy, axis=0) 
                basic_array = numpy.append(basic_array, mZ[1]+xz, axis=0) 
                known_case = True
            if (dipoleDirection.lower() == '011'):
                basic_array = numpy.append(Z[8]+Y[8]+lp, Z[8]+Y[8]+xy, axis=0)
                basic_array = numpy.append(basic_array, Z[8]+Y[8]+yz, axis=0)
                basic_array = numpy.append(basic_array, mZ[8]+mY[8]+xz, axis=0) 
                known_case = True

    if (not known_case):
        print "unhandled combination of lattice and dipole direction"
        print __doc__

    print basic_array

    bravais_lattice = []
    for i in range(latticeNumber):
        for j in range(latticeNumber):
            for k in range(latticeNumber):
                for l in range(len(basic_array)):
                    bravais_lattice.append(2*i*e1 + 2*j*e2 + 2*k*e3 + basic_array[l])

    outputFile = open(outputFileName, 'w')
    
    outputFile.write('<OpenMD version=2>\n')
    outputFile.write('  <MetaData>\n')
    outputFile.write('     molecule{\n')
    outputFile.write('       name = \"D\";\n')
    outputFile.write('       \n')
    outputFile.write('       atom[0]{\n')
    outputFile.write('         type = \"D\";\n')
    outputFile.write('         position(0.0, 0.0, 0.0);\n')
    outputFile.write('         orientation(0.0, 0.0, 0.0);\n')
    outputFile.write('       }\n')
    outputFile.write('     }\n')
    outputFile.write('     component{\n')
    outputFile.write('       type = \"D\";\n')
    outputFile.write('       nMol = '+ repr(len(bravais_lattice)) + ';\n')
    outputFile.write('     }\n')

    outputFile.write('     ensemble = NVE;\n')
    outputFile.write('     forceField = \"Multipole\";\n')

    outputFile.write('     cutoffMethod = \"shifted_force\";\n')
    outputFile.write('     electrostaticScreeningMethod = \"damped\";\n')

    outputFile.write('     cutoffRadius = 9.0;\n')
    outputFile.write('     dampingAlpha = 0.18;\n')
    outputFile.write('     statFileFormat = \"TIME|TOTAL_ENERGY|POTENTIAL_ENERGY|KINETIC_ENERGY|TEMPERATURE|PRESSURE|VOLUME|CONSERVED_QUANTITY|ELECTROSTATIC_POTENTIAL\";\n')
    outputFile.write('     dt = 1.0;\n')
    outputFile.write('     runTime = 1.0;\n')
    outputFile.write('     sampleTime = 1.0;\n')
    outputFile.write('     statusTime = 1.0;\n')
    outputFile.write('  </MetaData>\n')
    outputFile.write('  <Snapshot>\n')
    outputFile.write('    <FrameData>\n');
    outputFile.write("        Time: %.10g\n" % (0.0))
    
    Hxx = 2.0 * latticeConstant * latticeNumber
    Hyy = 2.0 * latticeConstant * latticeNumber
    Hzz = 2.0 * latticeConstant * latticeNumber
    
    outputFile.write('        Hmat: {{%d, 0, 0}, {0, %d, 0}, {0, 0, %d}}\n' % (Hxx, Hyy, Hzz))
    outputFile.write('    </FrameData>\n')
    outputFile.write('    <StuntDoubles>\n')
    sdFormat = 'pvqj'
    index = 0
    
    for i in range(len(bravais_lattice)):
        xcart = latticeConstant*(bravais_lattice[i][0])
        ycart = latticeConstant*(bravais_lattice[i][1])
        zcart = latticeConstant*(bravais_lattice[i][2])
        dx = bravais_lattice[i][3]
        dy = bravais_lattice[i][4]
        dz = bravais_lattice[i][5]

        uz = numpy.array([dx, dy, dz])
        uz = uz/numpy.linalg.norm(uz)
        
        uy = numpy.array([0.0, 1.0, 0.0])
        uy = uy - uz * numpy.vdot(uy, uz) / numpy.vdot(uz, uz) 
        uy = uy/numpy.linalg.norm(uy)

        ux = numpy.cross(uy, uz)

        RotMat = [ux, uy, uz]

        q = [0.0, 0.0, 0.0, 0.0]

        # RotMat to Quat code is out of OpenMD's SquareMatrix3.hpp code:

        t = RotMat[0][0] + RotMat[1][1] + RotMat[2][2] + 1.0
        
        if( t > 1e-6 ):
            s = 0.5 / math.sqrt( t )
            q[0] = 0.25 / s
            q[1] = (RotMat[1][2] - RotMat[2][1]) * s
            q[2] = (RotMat[2][0] - RotMat[0][2]) * s
            q[3] = (RotMat[0][1] - RotMat[1][0]) * s
        else:
            ad1 = RotMat[0][0]
            ad2 = RotMat[1][1]
            ad3 = RotMat[2][2]

            if( ad1 >= ad2 and ad1 >= ad3 ):
                s = 0.5 / math.sqrt( 1.0 + RotMat[0][0] - RotMat[1][1] - RotMat[2][2] )
                q[0] = (RotMat[1][2] - RotMat[2][1]) * s
                q[1] = 0.25 / s
                q[2] = (RotMat[0][1] + RotMat[1][0]) * s
                q[3] = (RotMat[0][2] + RotMat[2][0]) * s
            elif ( ad2 >= ad1 and ad2 >= ad3 ):
                s = 0.5 / math.sqrt( 1.0 + RotMat[1][1] - RotMat[0][0] - RotMat[2][2] )
                q[0] = (RotMat[2][0] - RotMat[0][2] ) * s
                q[1] = (RotMat[0][1] + RotMat[1][0]) * s
                q[2] = 0.25 / s
                q[3] = (RotMat[1][2] + RotMat[2][1]) * s
            else:
                s = 0.5 / math.sqrt( 1.0 + RotMat[2][2] - RotMat[0][0] - RotMat[1][1] )
                q[0] = (RotMat[0][1] - RotMat[1][0]) * s
                q[1] = (RotMat[0][2] + RotMat[2][0]) * s
                q[2] = (RotMat[1][2] + RotMat[2][1]) * s
                q[3] = 0.25 / s

        
        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))
        index = index+1

    outputFile.write("    </StuntDoubles>\n")
    outputFile.write("  </Snapshot>\n")
    outputFile.write("</OpenMD>\n")
    outputFile.close()

    outputFile.close()
    
def main(argv):
    
    arrayType = "A"
    haveOutputFileName = False
    latticeType = "fcc"
    dipoleDirection = "001"
    latticeNumber = 3
    latticeConstant = 4
    try:                                
        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="])  
    except getopt.GetoptError:           
        usage()                          
        sys.exit(2)                     
    for opt, arg in opts:                
        if opt in ("-h", "--help"):      
            usage()                     
            sys.exit()
        elif opt in ("-x", "--array-type-X"):
            arrayType = "X"
        elif opt in ("-y", "--array-type-Y"):
            arrayType = "Y"
        elif opt in ("-z", "--array-type-Z"):
            arrayType = "Z"
        elif opt in ("-b", "--array-type-B"):
            arrayType = "B"
        elif opt in ("-l", "--lattice"): 
            latticeType = arg
        elif opt in ("-d", "--direction"):
            dipoleDirection = arg
        elif opt in ("-c", "--constant"):
            latticeConstant = float(arg)
        elif opt in ("-n"):
            latticeNumber = int(arg)
        elif opt in ("-o", "--output-file"): 
            outputFileName = arg
            haveOutputFileName = True
    if (not haveOutputFileName):
        usage()
        print "No output file was specified"
        sys.exit()
        
    createLattice(latticeType, latticeNumber, latticeConstant, dipoleDirection, arrayType, outputFileName);

if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	1899
Committed:	Fri Jul 12 17:37:27 2013 UTC (12 years, 2 months ago) by gezelter
File size:	14021 byte(s)
Log Message:	Added the Z basic arrays
#	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 = (l1e1 + l2e2 + l3e3) + value e4
87			Yvec = (l1e1 + l2e2 + l3e3) + value e5
88			Zvec = (l1e1 + l2e2 + l3e3) + value e6
89			mXvec = (l1e1 + l2e2 + l3e3) - value e4
90			mYvec = (l1e1 + l2e2 + l3e3) - value e5
91			mZvec = (l1e1 + l2e2 + l3e3) - 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(2ie1 + 2je2 + 2ke3 + 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
300			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))
301			index = index+1
302
303			outputFile.write(" </StuntDoubles>\n")
304			outputFile.write(" </Snapshot>\n")
305			outputFile.write("</OpenMD>\n")
306	gezelter	1889	outputFile.close()
307
308	gezelter	1898	outputFile.close()
309
310	gezelter	1889	def main(argv):
311
312	gezelter	1899	arrayType = "A"
313	gezelter	1889	haveOutputFileName = False
314			latticeType = "fcc"
315			dipoleDirection = "001"
316			latticeNumber = 3
317			latticeConstant = 4
318			try:
319	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="])
320	gezelter	1889	except getopt.GetoptError:
321			usage()
322			sys.exit(2)
323			for opt, arg in opts:
324			if opt in ("-h", "--help"):
325			usage()
326			sys.exit()
327	gezelter	1899	elif opt in ("-x", "--array-type-X"):
328			arrayType = "X"
329			elif opt in ("-y", "--array-type-Y"):
330			arrayType = "Y"
331			elif opt in ("-z", "--array-type-Z"):
332			arrayType = "Z"
333	gezelter	1889	elif opt in ("-b", "--array-type-B"):
334	gezelter	1899	arrayType = "B"
335	gezelter	1889	elif opt in ("-l", "--lattice"):
336			latticeType = arg
337			elif opt in ("-d", "--direction"):
338			dipoleDirection = arg
339			elif opt in ("-c", "--constant"):
340			latticeConstant = float(arg)
341			elif opt in ("-n"):
342			latticeNumber = int(arg)
343			elif opt in ("-o", "--output-file"):
344			outputFileName = arg
345			haveOutputFileName = True
346			if (not haveOutputFileName):
347			usage()
348			print "No output file was specified"
349			sys.exit()
350
351	gezelter	1899	createLattice(latticeType, latticeNumber, latticeConstant, dipoleDirection, arrayType, outputFileName);
352	gezelter	1889
353			if __name__ == "__main__":
354			if len(sys.argv) == 1:
355			usage()
356			sys.exit()
357			main(sys.argv[1:])